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Abstract. Helioseismology studies the structure and dynamics of the Sun's in- 
terior by observing oscillations on the surface. These studies provide information 
about the physical processes that control the evolution and magnetic activity of 
the Sun. In recent years, helioseismology has made substantial progress towards 
the understanding of the physics of solar oscillations and the physical processes in- 
side the Sun, thanks to observational, theoretical and modeling efforts. In addition 
to the global seismology of the Sun based on measurements of global oscillation 
modes, a new field of local helioseismology, which studies oscillation travel times 
and local frequency shifts, has been developed. It is capable of providing 3D images 
of the subsurface structures and flows. The basic principles, recent advances and 
perspectives of global and local helioseismology are reviewed in this article. 



1 Introduction 

In 1926 in his book The Internal Constitution of the Stars Sir Arthur Stanley 
Eddington [1] wrote: 

At first sight it would seem that the deep interior of the sun and 
stars is less accessible to scientific investigation than any other region 
of the universe. Our telescopes may probe farther and farther into 
the depths of space; but how can we ever obtain certain knowledge 
of that which is hidden behind substantial barriers? What appliance 
can pierce through the outer layers of a star and test the conditions 
within? 

The answer to this question was provided a half a century later by helioseis- 
mology. Helioseismology studies the conditions inside the Sun by observing 
and analyzing oscillations and waves on the surface. The solar interior is 
not transparent to light but it is transparent to acoustic waves. Acoustic 
(sound) waves on the Sun are excited by turbulent convection below the vis- 
ible surface (photosphere) and travel through the interior with the speed of 
sound. Some of these waves are trapped inside the Sun and form resonant 
oscillation modes. The travel times of acoustic waves and frequencies of the 
oscillation modes depend on physical conditions of the internal layers (tem- 
perature, density, velocity of mass flows, etc). By measuring the travel times 
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and frequencies one can obtain information these condition. This is the ba- 
sic principle of helioseismology. Conceptually it is very similar to the Earth's 
seismology. The main difference is that the Earth's seismology studies mostly 
individual events, earthquakes, while helioseismology is based on the analysis 
of acoustic noise produced by solar convection. However, recently the local 
helioseismic techniques have been applied for ambient noise tomography of 
Earth's structures. The solar oscillations are observed in variations of inten- 
sity of solar images or, more commonly, in line-of-sight velocity of the surface 
elements, which is measured from the Doppler shift of spectral lines (Fig. 1). 
Variations caused by these oscillations are very small, much smaller than the 
noise produced by turbulent convection. Thus, their observation and analysis 
requires special procedures. 

Helioseismology is a relatively new discipline of solar physics and astro- 
physics. It has been developed over the past few decades by a large group 
of remarkable observers and theorists, and is continued being actively de- 
veloped. The history of helioseismology has been very fascinating, from the 
initial discovery of the solar 5-min oscillations and the initial attempts to un- 
derstand the physical nature and mechanism of these oscillations to detailed 
diagnostics of the deep interior and subsurface magnetic structures associ- 
ated with solar activity. This development was not straightforward. As this 
always happens in science controversial results and ideas provided inspiration 
for further more detailed studies. 

In a brief historical introduction, I describe some key contributions. It is 
very interesting to follow the line of discoveries that led to our current under- 
standing of the oscillations and helioseismology techniques. Then, I overview 
the basic concepts and results of helioseismology. The launch of the Solar 
Dynamics Observatory in 2010 opens a new era in helioseismology. The He- 
lioseismic and Magnetic Imager (HMI) instrument will provide uninterrupted 
high-resolution Doppler-shift and vector magnetogram data over the whole 
disk. These data will provide a complete information about the solar oscilla- 
tions and their interaction with solar magnetic fields. 

2 Brief history of helioseismology 

Solar oscillations were discovered in 1960 by Robert Leighton, Robert Noyers 
and George Simon [2] by analyzing series of Dopplergrams obtained at the 
Mt. Wilson Observatory. Instead of the expected turbulent behavior of the 
velocity field they found two distinct classes: large-scale horizontal cellular 
motions, which they called super granulation, and vertical quasi-periodic oscil- 
lations with a period of about 300 seconds (5 min) and a velocity amplitude 
of about 0.4 kms -1 . It turned out that these oscillations are the dominant 
vertical motion in the lower atmosphere (chromosphere) of the Sun. It is 
remarkable that they realized the diagnostic potential noting that these os- 
cillations "offer a new means of determining certain local properties of the 
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Fig. 1. a) Image of the line-of-sight (Doppler) velocity of the solar surface obtained 
by the Michelson Doppler Imager (MDI) instrument on board SOHO spacecraft on 
1997-06-19, 02:00 UT; b) Oscillations of the Doppler velocity, measured by MDI at 
the solar disk center in 12 CCD pixels separated by ~ 1.4 Mm on the Sun. 



solar atmosphere, such as the temperature, the vertical temperature gradient, 
or the mean molecular weight". They also pointed out that the oscillations 
might be excited in the Sun's granulation layer, and account for a part of the 
energy transfer from the convection zone into the chromosphere. 

This discovery was confirmed by other observers, and for several years 
it was believed that the oscillations represent transient atmospheric waves 
excited by granules, small convective cells on the solar surface, l-2x 10 3 km 
in size and 8—10 min lifetime. The physical nature of the oscillation at that 
time was unclear. In particular, the questions whether these oscillations are 
acoustic or gravity waves, and if they represent traveling or standing waves 
remained unanswered for almost a decade after the discovery. 

Pierre Mein [3] applied a two-dimensional Fourier analysis (in time and 
space) to observational data obtained by John Evans and his colleagues at 
the Sacramento Peak Observatory in 1962-65. His idea was to decompose 
the oscillation velocity field into normal modes. He calculated the oscillation 
power spectrum and investigated the relationship between the period and 
horizontal wavelength (or frequency-wavenumber diagram). From this analy- 
sis he concluded that the oscillations are acoustic waves that are stationary 
(evanescent) in the solar atmosphere. He also made a suggestion that the 
horizontal structure of the oscillations may be imposed by the convection 
zone below the surface. 

Mein's results were confirmed by Edward Frazier [4] who analyzed high- 
resolution spectrograms taken at the Kitt Peak National Observatory in 1965. 
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In the wavenumber-frequency diagram he noticed that in addition to the 
primary 5-min peak peak there is a secondary lower frequency peak, which 
was a new puzzle. 

This puzzle was solved by Roger Ulrich [5] who following the ideas of Mein 
and Frazier, calculated the spectrum of standing acoustic waves trapped in a 
layer below the photosphere. He found that these waves may exist only along 
discrete line in the wavenumber-frequency (k — w) diagram, and that the 
two peaks observed by Frazier correspond to the first two harmonics (nor- 
mal modes). He formulated the conditions for observing the discrete acoustic 
modes: observing runs must be longer longer one hour, must cover a suf- 
ficiently large region of, at least, 60,000 km in size; the Doppler velocity 
images must have a spatial resolution of 3,000 km, and be taken at least 
every 1 minute. 

At that time the observing runs were very short, typically, 30-40 min. Only 
in 1974-75 Franz-Ludwig Deubner [6] was able to obtain three 3-hour sets of 
observations using a magnetograph of the Fraunhofer Institute in Anacapri. 
He measured Doppler velocities along a ~ 220, 000 km line on the solar disk 
by scanning it periodically at 110 sec intervals with the scanning steps of 
about 700 km. The Fourier analysis of these data provided the frequency- 
wavenumber diagram with three or four mode ridges in the oscillation power 
spectrum that represents the squared amplitude of the Fourier components 
as a function of wavenumber and frequency. Deubner's results provided un- 
ambiguous confirmation of the idea that the 5-min oscillations observed on 
the solar surface represent the standing waves or resonant acoustic modes 
trapped below the surface. The lowest ridge in the diagram is easily identi- 
fied as the surface gravity wave because its frequencies depend only on the 
wavenumber and surface gravity. The ridge above is the first acoustic mode, 
a standing acoustic waves that have one node along the radius. The ridge 
above this corresponds to the second acoustic modes with two nodes, and so 
on. 

While these observations showed a remarkable qualitative agreement with 
Ulrich's theoretical prediction, the observed power ridges in the k — u dia- 
gram were systematically lower than the theoretical mode lines. Soon after, 
in 1975, Edward Rhodes, Ulrich and Simon [7] made independent observa- 
tions at the vacuum solar telescope at the Sacramento Peak Observatory 
and confirmed the observational results. They also calculated the theoretical 
mode frequencies for various solar models, and by comparing these with the 
observations determined the limits on the depth of the solar convection zone. 
This, probably, was the first helioseismic inference. 

However, it was believed that the acoustic (p) modes do not provide much 
information about the solar interior because detailed theoretical calculations 
of their properties by Hiroyashi Ando and Yoji Osaki [8] showed that while 
these mode are determined by interior resonances their amplitude (eigenfunc- 
tions) is predominantly concentrated close to the surface. Therefore, the main 
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focus was shifted to observations and analysis of global oscillations of the Sun 
with periods much longer than 5 min. This task was particularly important 
for explaining the observed deficit of high-energy solar neutrinos [9], which 
could be either due to a low temperature (or heavy element abundance - low 
metallicity) in the energy-generating core or neutrino oscillations. 

In 1975, Henry Hill, Tuck Stebbins and Tim Brown [10] reported on the 
detection of oscillations in their measurements of solar oblateness. The peri- 
ods of these oscillations were between 10 and 40 min. They suggested that 
the oscillation signals might correspond to global modes of the Sun. Inde- 
pendently, in 1976, two groups, led by Andrei Severny at the Crimean Ob- 
servatory [11] and George Isaak at the University of Birmingham [12] found 
long-period oscillations in global-Sun Doppler velocity signals. The oscillation 
with a period of 160 min was particularly prominent and stable. The ampli- 
tude of this oscillation was estimated close to 2 m/s. Later this oscillation 
was found in observations at the Wilcox Solar Observatory [13] and at the 
geographical South Pole [14]. Despite significant efforts to identify this oscil- 
lation among the solar resonant modes or to find a physical explanation these 
results remain a mystery. This oscillations lost the amplitude and coherence 
in the subsequent ground-based measurements and was not found in later 
observations from SOHO spacecraft [15]. The period of this oscillation was 
extremely close to 1/9 of a day, and likely was related to terrestrial observing 
conditions. 

Nevertheless, these studies played a very important role in development 
of helioseismology and emphasized the need for long-term stable and high- 
accuracy observations from the ground and space. Attempts to detect long- 
period oscillations (g-modes) still continue. However, the focus of helioseis- 
mology was shifted to accurate measurements and analysis of the acoustic 
p-modes discovered by Leighton. 

The next important step was made in 1979 by the Birmingham group [16]. 
They observed the Doppler velocity variations integrated over the whole Sun 
for about 300 hours (but typically 8 hours a day) at two observatories, Izana, 
on Tenerife, and Pic du Midi in the Pyrenees. In the power spectrum of 5-min 
oscillations they detected several equally space lines corresponding to global 
(low-degree) acoustic modes, radial, dipole and quadrupole. (In terms of the 
angular degree these are labeled as t = 0,1, and 2). Unlike, the previously 
observed local short horizontal wavelength acoustic modes these oscillations 
propagate into the deep interior and provide information about the structure 
of the solar core. The estimated frequency spacing between the modes was 
67.8 /iHz. This uniform spacing predicted theoretically by Yuri Vandakurov 
[17] in the framework of a general stellar oscillation theory corresponds to 
the inverse time that takes for acoustic waves to travel from the surface of 
the Sun through the center to the opposite side and come back. Thus, the 
frequency spacing immediately gives an important constraint on the internal 
structure of the Sun. A comparison with the solar models [18,19] showed 
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that the observed spectrum is consistent with the spectrum of solar models 
with low metallicity. This result was very exciting because if correct it would 
provide a solution to the solar neutrino problem. Thus, the determination of 
solar metallicity (or heavy element abundance) became a central problem of 
helioseismology. 

In the same year, 1974, Gerald Grec, Eric Fossat, and Martin Pomerantz 
[14] made 5-day continuous measurements at the Amundsen- Scott Station 
at the South Pole of the global oscillations and confirmed the Birmingham 
result. Also, they were able to resolve the fine structure of the oscillation 
spectrum and in addition to the main 67.8 /zHz spacing (large frequency 
separation) between the strongest peaks of £ = 1 and 2, observe a small 10-16 
/iHz splitting (small separation) between the £ = and 2, and £ = 1 and 
3 modes. The small separation is mostly sensitive to the central part of the 
Sun and provides additional diagnostic power. 

The comparison of the observed oscillation peaks in the frequency power 
spectra with the p-mode frequencies calculated for solar models showed that 
below the surface these oscillations correspond to the standing waves with a 
large number of nodes along the radius (or high radial order). The number 
of nodes is between 10 and 35, and it was difficult to determine the precise 
numbers for the observed modes. This created an uncertainty in the hclio- 
seismic determination of the heavy element abundance. Joergen Christensen- 
Dalsgaard and Douglas Gough [20] pointed out that while the South Pole 
and new Birmingham data favor solar models with normal metallicity the 
low metallicity models cannot be ruled out. 

The uncertainty was resolved three years later in 1983 when Tom Duvall 
and Jack Harvey [21] analyzed the Doppler velocity data measured with a 
photo-diode array in 200 positions along the North-South direction on the 
disk, and obtained the diagnostic k — oj diagram for acoustic modes of degree 
£, from 1 to 110. This allowed them to connect in the diagnostic diagram the 
global low-£ modes with the high-£ observed by Deubner. Since the correspon- 
dence of the ridges on Dcubncr's diagram to solar oscillation modes have been 
determined it was easy to identify the \ow-£ modes by simply counting the 
ridges corresponding to the low-^ frequencies. It turned out that the these 
modes are indeed in the best agreement with the normal metallicity solar 
model. This result had important implications for the solar neutrino prob- 
lem because it strongly indicated that the observed deficit of solar neutrinos 
was not due to a low abundance of heavy elements on the Sun but because 
of changes in neutrino properties (neutrino oscillations) on their way from 
the energy-generating core to the Earth. This was later confirmed by direct 
measurements of solar neutrino properties [22] . 

It was also important that the definite identification of the observed solar 
oscillations in terms of normal oscillation modes provided a solid founda- 
tion for developing diagnostic methods of helioseismology based on the well- 
developed mathematical theory of non-radial oscillations of stars [23-25]. 
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This theory provided means for calculating eigenfrequencies and eigenfunc- 
tions of normal modes for spherically symmetric stellar models. Mathemati- 
cally, the problem is reduced to solving a non-linear eigenvalue problem for a 
fourth-order system of differential equations. This system has two sequences 
of eigenvalues corresponding to p- and g-modes, and also a degenerate solu- 
tion, corresponding to f-modes (surface gravity waves). The effects of rota- 
tion, asphericity and magnetic fields are usually small and considered by a 
perturbation theory [26-29] . 

An important prediction of the oscillation theory is that rotation causes 
splitting of normal mode frequencies. Without rotation, the normal mode fre- 
quencies are degenerate with respect to the azimuthal wavenumbcr, m, that 
is the modes of the angular degree, /, and radial order, n, have the same fre- 
quencies irrespective of the azimuthal (longitudinal) wavelength. The stellar 
rotation removes this degeneracy. Obviously, it does not affect the axisym- 
metrical (m=0) modes, but the frequencies of non-axisymmetrical modes are 
split. Generally, these modes can be represented as a superposition of two 
waves running around a star in two opposite directions (prograde and retro- 
grade waves) . Without rotation, these modes have the same frequencies and, 
thus, the same phase speed. In this case, they form a standing wave. However, 
rotation increases the speed of the prograde wave and decreases the speed 
of retrograde wave. This results in an increase of the eigenfrequency of the 
prograde mode, and a frequency decrease of the retrograde mode. This phe- 
nomenon is similar to frequency shifts due to the Doppler effect. It is called 
rotational frequency splitting. 

The rotational frequency splitting was first observed by Ed Rhodes, Roger 
Ulrich and Franz Deubner [30-32]. These measurements provided first evi- 
dence that the rotation rate of the Sun is not uniform but increases with 
depth. The rotational splitting was initially measured for high-degree modes, 
but then the measurements were extended to medium- and low-degree range 
by Tom Duvall and Jack Harvey [33,34], who made a long continuous series 
of helioseismology observations at the South Pole. The internal differential 
rotation law was determined from the data of Tim Brown and Cherilynn Mor- 
row [35]. It was found that the differential latitudinal rotation is confined in 
the convection zone, and that the radiative interior rotates almost uniformly, 
and also slower in the equatorial region than the convective envelope [36,37]. 
Such rotation law was not expected from theories of stellar rotation, which 
predicted that the stellar cores rotate faster than the envelopes [38]. The 
knowledge of the Sun's internal rotation law is of particular importance for 
understanding the dynamo mechanism of magnetic field generation [39]. 

It became clear that for long uninterrupted observations are essential for 
accurate inferences of the internal structure and rotation of the Sun. There- 
fore, the observational programs focused on development of global helioseis- 
mology networks, GONG [40] and BiSON [41,42], and also the Solar and He- 
liospheric Observatory (SOHO) space mission [43]. These projects provided 
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almost continuous coverage for helioseismic observations and also stimulated 
development of new sophisticated data analysis and inversion techniques. 

In addition, the Michelson Doppler Imager (MDI) instrument on SOHO 
[44] and the GONG+ network upgraded to higher spatial resolution [45] 
provided excellent opportunities for developing local helio seismology, which 
provides tools for three-dimensional imaging of the solar interior. The lo- 
cal helioseismology methods are based on measurements of local oscillation 
properties, such as frequency shift in local areas or variations of travel times. 

The idea of using the local frequency shifts for inferring the subsurface 
flows was suggested by Douglas Gough and Juri Toomre in 1983 [46]. The 
method is now called ring-diagram analysis [47], because the dispersion re- 
lation of solar oscillations forms rings in horizontal wavenumber plane at a 
given frequency. It measures shifts of these rings, which are then converted 
into frequency shifts. Ten years later, Tom Duvall and his colleagues [48] 
introduced time- distance helioseismology method. In this method, they sug- 
gested to measure travel times of acoustic waves from a cross-covariance 
function of solar oscillations. This function is obtained by cross-correlating 
oscillation signals observed at two different points on the solar surface for 
various time lags. When the time lag in the calculations coincides with the 
travel time of acoustic waves between these points the cross-covariance func- 
tion shows a maximum. This method provided means for developing acoustic 
tomography techniques [49,50] for imaging 3D structures and flows with the 
high-resolution comparable to the oscillation wavelength. These and other 
methods of local area helioseismology [51,52] have provided important re- 
sults on the convective and large-scale flows, and also on the structure and 
evolution of sunspots and active regions. Their development continues. 

The SOHO mission and the GONG network were primarily designed for 
observing solar oscillation modes of low- and medium-degree, needed for 
global helioseismology. Local helioseismology requires high-resolution obser- 
vations of high-degree modes. Because of the telemetry constraints such data 
are available uninterruptedly from the MDI instrument on SOHO only for 
2 months every year. These data provided only snapshots of the subsurface 
structures and dynamics associated with the solar activity. In order to fully 
investigate the evolving magnetic activity of the Sun, a new space mission 
Solar Dynamics Observatory (SDO) was launched on February 11, 2010. It 
carries Helioseismic and Magnetic Imager (HMI) instrument, which will pro- 
vide continuous 4096x4096-pixcl full-disk images of solar oscillations. These 
data will open new opportunities for investigation the solar interior by local 
helioseismology [53]. 

In the modern helioseismology, a very important role is played by nu- 
merical simulations. Both, global and local helioseismology analysis employ 
relatively simple for fitting the observational data and performing inversions 
of the fitted frequencies and travel times. For instance, the global helioseis- 
mology methods assume that the structures and flows on the Sun are axisym- 



Advances in Global and Local Helioseismology 9 



metrical and infer only the axisymmetrical components of the sound speed 
and velocity field. The local helioseismology methods are based on a simplified 
physics of wave propagation on the Sun. The ring-diagram analysis makes an 
assumption that that the perturbations and flows are horizontally uniform 
within the area used for calculating the wave dispersion relation, 5-15 hclio- 
graphic degrees, while a typical size of sunspots is about 1-2 degrees. Most of 
the time-distance helioseismology inversions are based on a ray-path approx- 
imation and ignore the finite wavelength effects that become important at 
small scales, comparable with the wavelength. Also, all the methods, global 
and local, do not take into account many effects of solar magnetic fields. Prop- 
erties of solar oscillations dramatically change in regions of strong magnetic 
field. In particular, the excitation of oscillations is suppressed in sunspots 
because the strong magnetic field inhibits convection that drives the oscilla- 
tions. The magnetic stresses may cause anisotropy of wave speed and lead 
to transformation of acoustic waves into various MHD type waves. These 
and other effects have to be investigated and taken into account in the data 
analysis and inversion procedures. Because of the complexity, these processes 
can be fully investigated only numerically. The numerical simulations of sub- 
surface solar convection and oscillations were pioneered by Robert Stein and 
Ake Nordlund [54] . These 3D radiative MHD simulations include all essential 
physics and provide important insights into the physical processes below the 
visible surface and also artificial data for helioseismology testing. This type 
of so-called "realistic" simulations has been used for testing time-distance 
helioseismology inferences [55] , and continues being developed using modern 
turbulence models [56] . In addition, for testing various aspects of wave propa- 
gation and interaction with magnetic fields are studied by solving numerically 
linearized MHD equations (e.g. [57-59]). The numerical simulations become 
an important tool for verification and testing of the helioseismology methods 
and inferences. 

3 Basic properties of solar oscillations 
3.1 Oscillation power spectrum 

The theoretical spectrum of solar oscillation modes shown in Fig. 2 covers 
a wide range of frequencies and angular degrees. It includes oscillations of 
three types: acoustic (p) modes, surface gravity (f) modes and internal grav- 
ity (g) modes. In this spectrum, the modes are organized a series of curves 
corresponding to different overtones of non-radial modes, which are charac- 
terized by the number of nodes along the radius (or by the radial order, n). 
The angular degree, I, of the corresponding spherical harmonics describes 
the horizontal wave number (or inverse horizontal wavelength) . The p-modes 
cover the frequency range from 0.3 to 5 mHz (or from 3 to 55 min in oscillation 
periods). The low frequency limit corresponds to the first radial harmonic, 
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Fig. 2. Theoretical frequencies of solar oscillation modes calculated for a standard 
solar model for a range of angular degree I from to 100, and for the frequency 
range from 0.2 mHz to 5 mHz. The solid curves connect modes corresponding to 
the different oscillation overtones (radial orders) The dashed grey horizontal line 
indicate the low-frequency observational limit: only the modes above this line have 
been reliably observed. The right panel shows an area of the avoided crossing of f- 
and g-modes (indicated by the gray dashed circle in left panel). 



and the upper limit is set by the acoustic cut-off frequency of the solar at- 
mosphere. The g-modes frequencies have an upper limit corresponding to the 
maximum Brunt- Valsala frequency (~ 0.45 mHz) in the radiative zone and 
occupy the low-frequency part of the spectrum. The intermediate frequency 
range of 0.3-0.4 mHz at low angular degrees is a region of mixed modes. These 
modes behave like g-modes in the deep interior and like p-modes in the outer 
region. The apparent crossings in this diagram are not the actual crossings: 
the mode branches become close in frequencies but do not cross each other. 
At these points the mode exchange their properties, and the mode branches 
are diverted. For instance, the f-mode ridge stays above the g-mode lines. A 
similar phenomenon is known in quantum mechanics as avoided crossing. 

So far, only the upper part of the solar oscillation spectrum is observed. 
The lowest frequencies of detected p- and f- modes are of about 1 mHz. At 
lower frequencies the mode amplitudes decrease below the noise level, and 
become unobservable. There have been several attempts to identify low- 
frequency p-modes or even g-modes in the noisy spectrum, but so far these 
results are not convincing. 

The observed power spectrum is shown in Fig. 3. The lowest ridge is the 
f-mode, and the other ridges are p-modes of the radial order, n, starting 
from n=l. The ridges of the oscillation modes disappear in the convective 
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Fig. 3. Power spectrum obtained from a 6-day long time series of solar oscillation 
data from the MDI instrument on SOHO in 1996 (y is the cyclic frequency of the 
oscillations, I is the angular degree, \h is the horizontal wavelength in megameters). 



noise at frequencies below 1 mHz. The power spectrum is obtained from the 
SOHO/MDI data, representing 1024xl024-pixel images of the line-of-sight 
(Doppler) velocity of the solar surface taken every minute without interrup- 
tion. When the oscillations are observed in the integrated solar light ("Sun- 
as-a-star") then only the modes of low angular degree are detected in the 
power spectrum (Fig. 4). These modes have a mean period of about 5 min, 
and represent p-modes of high radial order n modes. The n-values of these 
modes can be determined by tracing in Fig. 3 the the high-n ridges of the 
high-degree modes into the low-degree region. This provides unambiguous 
identification of the low-degree solar modes. Obviously, the mode identifica- 
tion is much more difficult for spatially unresolved oscillations of other stars. 

3.2 Excitation by turbulent convection 

Observations and numerical simulations have shown that solar oscillations 
are driven by turbulent convection in a shallow subsurface layer with a su- 
peradiabatic stratification, where convective velocities are the highest. How- 
ever, details of the stochastic excitation mechanism are not fully established. 
Solar convection in the superadiabatic layer forms small-scale granulation 
cells. Analysis of the observations and numerical simulations has shown that 
sources of solar oscillations are associated with strong downdrafts in dark in- 
tergranular lanes [60]. These downdrafts are driven by radiative cooling and 
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Fig. 4. Power spectral density (PSD) of low-degree solar oscillations, obtained 
from the integrated light observations (Sun-as-a-star) by the GOLF instrument 
on SOHO, from 11/04/1996 to 08/07/2008. 



may reach near-sonic velocity of several km/s. This process has features of 
convective collapse [61]. 

Calculations of the work integral for acoustic modes using the realistic 
numerical simulations of Stein and Nordlund [62] have shown that the prin- 
cipal contribution to the mode excitation is provided by turbulent Reynolds 
stresses and that a smaller contribution comes from non-adiabatic pressure 
fluctuations. Because of the very high Reynolds number of the solar dynamics 
the numerical modeling requires an accurate description of turbulent dissi- 
pation and transport on the numerical subgrid scale. The recent radiative 
hydrodynamics modeling using the Large-Eddy Simulations (LES) approach 
and various subgrid scale (SGS) formulations [56] showed that among these 
formulations the most accurate description in terms of the reproducing the 
total amount of the stochastic energy input to the acoustic oscillations is 
provided by a dynamic Smagorinsky model [63,64] (Fig. 5a). 

As we have pointed out, the observations show that the modal lines in 
the oscillation power spectrum are not Lorentzian but display a strong asym- 
metry [67,68]. Curiously, the asymmetry has the opposite sense in the power 
spectra calculated from Dopplcr velocity and intensity oscillations. The asym- 
metry itself can be easily explained by interference of waves emanated by a 
localized source [69], but the asymmetry reversal is surprising and indicates 
complicated radiative dynamics of the excitation process. The reversal has 
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Fig. 5. a) Comparison of observed and calculated rate of stochastic energy in- 
put to modes for the entire solar surface (erg s _1 ). Different curves show the 
numerical simulation results obtained for 4 turbulence models: hyperviscosity 
(solid), enhanced hyperviscosity (dots), Smagorinsky (dash-dots), and dynamic 
model (dashes). Observed distributions: circles SoHo-GOLF, squares BISON, and 
triangles GONG for I — 1 [65]. b) Logarithm of the work integrand in units of 
erg cm~ 2 s _1 ), as a function of depth and frequency for numerical simulations 
with the dynamic turbulence model [66]. 



been attributed to a correlated noise contribution to the observed intensity 
oscillations [70], but the physics of this effect is still not fully understood. 
However, it is clear that the line shape of the oscillation modes and the phase- 
amplitude relations of the velocity and intensity oscillations carry substantial 
information about the excitation mechanism and, thus, require careful data 
analysis and modeling. 

3.3 Line asymmetry and pseudo-modes 

Figure 6 shows the power spectrum for oscillations of the angular degree, 
I = 200, obtained from the SOHO/MDI Doppler velocity and intensity data 
[70]. The line asymmetry is apparent, particularly, at low frequencies. In the 
velocity spectrum, there is more power in the low-frequency wings than in 
the high-frequency wings of the spectral lines. In the intensity spectrum, the 
distribution of power is reversed. The data also show that the asymmetry 
varies with frequency. It is the strongest for the f-mode and low-frequency p- 
mode peaks. At higher frequencies the peaks become more symmetrical, and 
extend well above the acoustic cut-off frequency (Eq. 51), which is ~ 5 — 5.5 
mHz. 

Acoustic waves with frequencies below the cut-off frequency are com- 
pletely reflected by the surface layers because of the steep density gradient. 
These waves are trapped in the interior, and their frequencies are deter- 
mined by the resonant conditions, which depend on the solar structure. But 
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Velocity and intensity spectra from SOHO/MD] 
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Fig. 6. Power spectra of / = 200 modes obtained from SOHO/MDI observations of 
a) Doppler velocity, b) continuum intensity [70] . 



the waves with frequencies above the cut-off frequency escape into the so- 
lar atmosphere. Above this frequency the power spectrum peaks correspond 
to so-called "pseudo- modes". These are caused by constructive interference 
of acoustic waves excited by the sources located in the granulation layer and 
traveling upward, and by the waves traveling downward, reflected in the deep 
interior and arriving back to the surface. Frequencies of these modes are no 
longer determined by the resonant conditions of the solar structure. They 
depend on the location and properties of the excitation source (" source reso- 
nance"). The pseudo-mode peaks in the velocity and intensity power spectra 
arc shifted relative to each other by almost a half- width. They are also slightly 
shifted relative to the normal mode peaks although they look like a contin- 
uation of the normal-mode ridges in Figs lb and 4a. This happens because 
the excitation sources are located in a shallow subsurface layer, which is very 
close to the reflection layers of the normal modes. Changes in the frequency 
distributions below and above the acoustic cut-off frequency can be easily 
noticed by plotting the frequency differences along the modal ridges. 

The asymmetrical profiles of normal-mode peaks are also caused by the 
localized excitation sources. The interference signal between acoustic waves 
traveling from the source upwards and the waves traveling from the source 
downward and coming back to the surface after the internal reflection de- 
pends on the wave frequency. Depending on the multipole type of the source 
the interference signal can be stronger at frequencies lower or higher than the 
resonant normal frequencies, thus resulting in asymmetry in the power dis- 
tribution around the resonant peak. Calculations of Nigam ct al. [70] showed 
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Fig. 7. a) The oscillation power spectrum from Hinode Call H line observations, 
b) The phase shift between Call H and G-band (units are in radians) [71]. 



that the asymmetry observed in the velocity spectra and the distribution of 
the pseudo-mode peaks can be explained by a composite source consisting 
of a monopole term (mass term) and a dipole term (force due to Reynolds 
stress) located in the zone of superadiabatic convection at a depth of ~ 100 
km below the photosphere. In this model, the reversed asymmetry in the in- 
tensity power spectra is explained by effects of a correlated noise added to the 
oscillation signal through fluctuations of solar radiation during the excitation 
process. Indeed, if the excitation mechanism is associated with the high-speed 
turbulent downdrafts in dark lanes of granulation the local darkening con- 
tributes to the intensity fluctuations caused by excited waves. The model 
also explains the shifts of pseudo-mode frequency peaks and their higher am- 
plitude in the intensity spectra. The difference between the correlated and 
uncorrelated noise is that the correlated noise has some phase coherence with 
the oscillation signal, while the uncorrelated noise has no coherence. 

While this scenario looks plausible and qualitatively explains the main 
properties of the power spectra details of the physical processes are still un- 
certain. In particular, it is unclear whether the correlated noise affects only 
the intensity signal or both the intensity and velocity. It has been suggested 
that the velocity signal may have a correlated contribution due to convec- 
tive overshoot [72]. Attempts to estimate the correlated noise components 
from the observed spectra have not provided conclusive results [73,74]. Re- 
alistic numerical simulations [75] have reproduced the observed asymmetries 
and provided an indication that radiation transfer plays a critical role in the 
asymmetry reversal. 

Recent high-resolution observations of solar oscillation simultaneously in 
two intensity filters, in molecular G-band and Call H line, from the Hinode 
space mission [76,77] revealed significant shifts in frequencies of pseudo-modes 
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observed in the Call H and G-band intensity oscillations [71] . The phase of the 
cross-spectrum of these oscillations shows peaks associated with the p-modc 
lines but no phase shift for the f-mode (Fig. 7b). The p-mode properties can 
be qualitatively reproduced in a simple model with a correlated background 
if the correlated noise level in the Ca II H data is higher than in the G-band 
data [71]. Perhaps, the same effect can explain also the frequency shift of 
pseudo-modes. The Call H line is formed in the lower chromosphere while 
the G-band signal comes from the photosphere. But how this may lead to 
different levels of the correlated noise is unclear. 

The Hinode results suggest that multi-wavelength observations of solar 
oscillations, in combination with the traditional intensity-velocity observa- 
tions, may help to measure the level of the correlated background noise and 
to determine the type of wave excitation sources on the Sun. This is impor- 
tant for understanding the physical mechanism of the line asymmetry and 
for developing more accurate models and fitting formulae for determining the 
mode frequencies [78]. 

In addition, Hinode provided observations of non-radial acoustic and sur- 
face gravity modes of very high angular degree. These observations show that 
the oscillation ridges are extended up to I ~ 4000 (Fig. 7a). In the high-degree 
range, I > 2500 frequencies of all oscillations exceed the acoustic cut-off fre- 
quency. The line width of these oscillations dramatically increases, probably 
due to strong scattering on turbulence [79,80]. Nevertheless, the ridge struc- 
ture extending up to 8 mHz (Nyquist frequency of these observations) is 
quite clear. Although the ridge slope clearly changes at the transition from 
the normal modes to the pseudo-modes. 

3.4 Magnetic effects: sunspot oscillations and acoustic halos 

In general, the main factors causing variations in oscillation properties in 
magnetic regions, can be divided in two types: direct and indirect. The direct 
effects are due to additional magnetic restoring forces that can change the 
wave speed and may transform acoustic waves into different types of MHD 
waves. The indirect effects are caused by changes in convective and ther- 
modynamic properties in magnetic regions. These include depth-dependent 
variations of temperature and density, large-scale flows, and changes in wave 
source distribution and strength. Both direct and indirect effects may be 
present in observed properties such as oscillation frequencies and travel times, 
and often cannot be easily disentangled by data analyses, causing confusions 
and misinterpretations. Also, one should keep in mind that simple models of 
MHD waves derived for various uniform magnetic configurations and without 
stratification or with a polytropic stratification may not provide correct ex- 
planations to solar phenomena. In this situation, numerical simulations play 
an important role in investigations of magnetic effects. 

Observed changes of oscillation amplitude and frequencies in magnetic 
regions are often explained as a result of wave scattering and conversion into 



Advances in Global and Local Helioseismology 



17 



various MHD modes. However, recent numerical simulations helped us to 
understand that magnetic fields not only affect the wave dispersion properties 
but also the excitation mechanism. In fact, changes in excitation properties 
of turbulent convection in magnetic regions may play a dominant role in 
observed phenomena. 

Sunspot oscillations For instance, it is well-known that the amplitude of 5- 
min oscillations is substantially reduced in sunspots. Observations show that 
more waves are coming into the sunspot than going out of the sunspot area 
(e.g. [81]). This is often attributed to absorption of acoustic waves in magnetic 
field due to conversion into slow MHD modes traveling along the field lines 
(e.g. [82]). However, since convective motions are inhibited by the strong mag- 
netic field of sunspots, the excitation mechanism is also suppressed. Three- 
dimensional numerical simulations of this effect have shown that the reduc- 
tion of acoustic emissivity can explain at least 50% of the observed power 
deficit in sunspots (Fig. 8) [83]. 




Fig. 8. a) Line-of-sight magnetic field map of a sunspot (AR8243); b) oscillation 
amplitude map; c) profiles of rms oscillation velocities at frequency 3.65 mHz for 
observations (thick solid curves) and simulations (dashed curves); the thin solid 
curve shows the distribution of the simulated source strength [83]. 



Another significant contribution comes from the amplitude changes caused 
by variations in the background conditions. Inhomogeneities in the sound 
speed may increase or decrease the amplitude of acoustic wave traveling 
through these inhomogeneities. Numerical simulations of MHD waves us- 
ing magnetostatic sunspot models show that the amplitude of acoustic waves 
traveling through sunspot decreases when the wave is inside sunspot and then 
increases when the wave comes out of sunspot [84] . Simulations with multiple 
random sources show that these changes in the wave amplitude together with 
the suppression of acoustic sources can explain the whole observed deficit of 
the power of 5-min oscillations. Thus, the role of the MHD mode conver- 
sion may be insignificant for explaining the power deficit of 5-min photo- 
spheric oscillations in sunspots. However, the mode conversion is expected to 



18 Alexander G. Kosovichev 



be significant higher in the solar atmosphere where magnetic forces become 
dominant. 

We should note that while the 5-min oscillations in sunspots come mostly 
from outside sources there are also 3-min oscillations, which are probably 
intrinsic oscillations of sunspots. The origin of these oscillations is not yet 
understood. They are probably excited by a different mechanism operating 
in strong magnetic field. 




Fig. 9. Call H intensity image from Hinode observations (top-left) and the cor- 
responding power maps from Call H intensity data in five frequency intervals of 
active region NOAA 10935. The field of view is 100 arcsec square in all the panels. 
The power is displayed in logarithmic greyscaling [85]. 



Hinode observations added new puzzles to sunspot oscillations. Figure 3.4 
shows a sample Ca II H intensity and the relative intensity power maps aver- 
aged over 1 mHz intervals in the range from 1 mHz to 7 mHz with logarithmic 
greyscaling [85]. In the CallH power maps, in all the frequency ranges, there 
is a small area (~ 6 arcsec in diameter) near the center of the umbra where 
the power was suppressed. This type of 'node' has not been reported be- 
fore. Possibly, the stable high-resolution observation made by Hinode/SOT 
was required to find such a tiny node, although analysis of other sunspots 
indicates that probably only a particular type of sunspots, e.g., round ones 
with axisymmetric geometry, exhibit such node-like structure. Above 4 mHz 
in the CallH power maps, power in the umbra is remarkably high. In the 
power maps averaged over narrower frequency range (0.05 mHz wide, not 
shown), the region with high power in the umbra seems to be more patchy. 
This may correspond to elements of umbral flashes, probably caused by over- 
shooting convective elements [86]. The CallH power maps show a bright 
ring in the penumbra at lower frequencies. It probably corresponds to the 
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running penumbral waves. The power spectrum in the umbra has two peaks: 
one around 3 mHz and the other around 5.5 mHz. The high-frequency peak 
is caused by the oscillations that excited only in the strong magnetic field of 
sunspots. The origin of these oscillations is not known yet. 



A^7^^ne^j^iGL^aoneti^Kl^ OsctEtatkm Power Mirp. v 2.5-3.B mHz OiiCtLLuttou Power Map, v 5.3-6,4 mHz 




Fig. 10. a) Line-of-sight magnetic field map of active region NOAA 9787 observed 
from SOHO/MDI on Jan. 24, 2002 and averaged over a 3- hour period; b) oscilla- 
tion power map from Doppler velocity measurements for the same period in the 
frequency 2.5-3.8 mHz; c) power map for 5.3-6.4 mHz. 



Acoustic halos In moderate field regions, such as plages around sunspot 
regions, observations reveal enhanced emission at high frequencies, 5-7 mHz, 
(with period <~ 3 min) [87]. Sometimes this emission is called the "acous- 
tic halo" (Fig. 10c). There have been several attempts to explain this effect 
as a result of wave transformation or scattering in magnetic structures (e.g. 
[88,89]). However, numerical simulations show that magnetic field can change 
the excitation properties of solar granulation resulting in an enhanced high- 
frequency emission. In particular, the radiative MHD simulations of solar 
convection [66] in the presence of vertical magnetic field have shown that 
the magnetic field significantly changes the structure and dynamics of gran- 
ulations, and thus the conditions of wave excitation. In magnetic field the 
granules become smaller, and the turbulence spectrum is shifted towards 
higher frequencies. This is illustrated in Figure 11, which shows the frequency 
spectrum of the horizontally averaged vertical velocity. Without a magnetic 
field the turbulence spectrum declines sharply at frequencies above 5 mHz, 
but in the presence of magnetic field it develops a plateau. In the plateau 
region characteristic peaks (corresponding to the "pseudo-modes") appear 
in the spectrum for moderate magnetic field strength of about 300-600 G. 
These peaks may explain the effect of the " acoustic halo" . Of course, more 
detailed theoretical and observational studies are required to confirm this 
mechanism. In particular, multi-wavelength observations of solar oscillations 
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at several different heights would be important. Investigations of the excita- 
tion mechanism in magnetic regions is also important for interpretation of 
the variations of the frequency spectrum of low-degree modes on the Sun, 
and for asteroseismic diagnostics of stellar activity. 




Fig. 11. Power spectra of the horizontally averaged vertical velocity at the visible 
surface for different initial vertical magnetic fields. The peaks on the top of the 
smooth background spectrum of turbulent convection represent oscillation modes: 
the sharp asymmetric peaks below 6 mHz are resonant normal modes, while the 
broader peaks above 6 mHz, which become stronger in magnetic regions, correspond 
to pseudo-modes. [66] 



3.5 Impulsive excitation: sunquakes 

"Sunquakes", the hclioseismic response to solar flares, are caused by strong 
localized hydrodynamic impacts in the photosphere during the flare impul- 
sive phase. The helioseismic waves have been observed directly as expanding 
circular-shaped ripples in SOHO/MDI Dopplcrgrams [90] (Fig. 12). 

These waves can be detected in Dopplcrgram movies and as a charac- 
teristic ridge in time-distance diagrams (Fig. 13a), [90-93], or indirectly by 
calculating integrated acoustic emission [94-96]. Solar flares are sources of 
high-temperature plasma and strong hydrodynamic motions in the solar at- 
mosphere. Perhaps, in all flares such perturbations generate acoustic waves 
traveling through the interior. However, only in some flares is the impact 
sufficiently localized and strong to produce the seismic waves with the am- 
plitude above the convection noise level. It has been established in the initial 
July 9, 1996, flare observations [90] that the hydrodynamic impact follows 
the hard X-ray flux impulse, and hence, the impact of high-energy electrons. 
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Fig. 12. Observations of the seismic response ( "sunquakes" ) of the solar flare 
of 9 July, 1996, showing a sequence of Doppler-velocity images, taken by the 
SOHO/MDI instrument. The signal of expanding ripples is enhanced by a factor 4 
in the these images. 




Fig. 13. a) The time-distance diagram of the seismic response to the solar flare of 9 
July, 1996. b) Illustration of acoustic ray paths of the flare-excited waves traveling 
through the Sun. 



A characteristic feature of the seismic response in this flare and several 
others [91-93] is anisotropy of the wave front: the observed wave amplitude is 
much stronger in one direction than in the others. In particular, the seismic 
waves excited during the October 28, 2003, 16 July, 2004, flare of 15 January, 
2005 flare had the greatest amplitude in the direction of the expanding flare 
ribbons (Fig. 14). The wave anisotropy can be attributed to the moving source 
of the hydrodynamic impact, which is located in the flare ribbons [91,93,97]. 
The motion of flare ribbons is often interpreted as a result of the magnetic 
reconnection processes in the corona. When the reconnection region moves up 
it involves higher magnetic loops, the footpoints of which are further apart. 
The motion of the footpoints of impact of the high-energy particles is par- 
ticularly well observed in the SOHO/MDI magnetograms showing magnetic 
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Fig. 14. Observations of the seismic response of the Sun ( "sunquakes" ) to two 
solar flares: a-c) X3 of 16 July, 2004, and d-f) XI flare of 15 January, 2005. The 
left panels show a superposition of MDI white-light images of the active regions 
and locations of the sources of the seismic waves determined from MDI Doppler- 
grams, the middle column shows the seismic waves, and the right panels show the 
time-distance diagrams of these events. The thin yellow curves in the right panels 
represent a theoretical time-distance relation for helioseismic waves for a standard 
solar model. [93] 



transients moving with supersonic speed, in some cases [92]. Of course, there 
might be other reasons for the anisotropy of the wave front, such as inho- 
mogeneities in temperature, magnetic field, and plasma flows. However, the 
source motion seems to be a key factor. 

Therefore, we conclude that the seismic wave was generated not by a 
single impulse but by a series of impulses, which produce the hydrodynamic 
source moving on the solar surface with a supersonic speed. The seismic effect 
of the moving source can be easily calculated by convolving the wave Green's 
function with a moving source function. The results of these calculations a 
strong anisotropic wavefront, qualitatively similar to the observations [97]. 
Curiously, this effect is quite similar to the anisotropy of seismic waves on 
Earth, when the earthquake rupture moves along the fault. Thus, taking into 
account the effects of multiple impulses of accelerated electrons and moving 
source is very important for sunquake theories. The impulsive sunquake oscil- 
lations provide unique information about interaction of acoustic waves with 
sunspots. Thus, these effects must be studied in more detail. 
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4 Global helioseismology 
4.1 Basic equations 

A simple theoretical model of solar oscillations can be derived using the 
following assumptions: 

1. linearity: v/c « 1, where v is velocity of oscillating elements, c is the 
speed of sound; 

2. adiabaticity: dS/dt — 0, where S is the specific entropy; 

3. spherical symmetry of the background state; 

4. magnetic forces and Reynolds stresses are negligible. 

The basic governing equations are derived from the conservation of mass, 
momentum, energy and the Newton's gravity law. The conservation of mass 
(continuity equation) assumes that the rate of mass change in a fluid element 
of volume V is equal to the mass flux through the surface of this element (of 
area A): 

^- t J v PdV = -J^pvda = -J^V(pv)dV, (1) 
where p is the mass density. Then, 

^+V(H = 0, (2) 
or in terms of the material derivative dp/dt = dp/dt + v ■ Vp: 

§ + ,V„ = 0. (3) 

The momentum equation (conservation of momentum of a fluid element) is: 

dv _ , „ . 

p^ = -VP + pg, (4) 

where P is pressure, g is the gravity acceleration, which can be expressed 
in terms of gravitational potential <P: g — V<£, dv/dt = dv/dt + v ■ Vv 
is the material derivative for the velocity vector. The adiabaticity equation 
(conservation of energy) for a fluid element is: 

or 

dP 2 dp 

-t= C Tt' (6) 

where c 2 = "fP/p is the squared adiabatic sound speed. The gravitational 
potential is calculated from the Poisson equation: 

W 2 <P = 4nGp. (7) 
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Now, we consider small perturbations of a stationary spherically symmet- 
rical star in hydrostatic equilibrium: 

v = 0, p = P o(r), P = P (r). 

If £(t) is a vector of displacement of a fluid element then velocity v of this 
clement: 

dt dt w 
Perturbations of scalar variables, p, P, can be of two general types: Eulerian 
(denoted with prime symbol), at a fixed position r: 

p(r,t) = po(r) +p'(r,t), 

and Lagrangian, measured in the moving element (denoted with 6): 

6p(r + £)=p (r) + 6p(r,t). (9) 

The Eulerian and Lagrangian perturbations are related to each other: 

5p = p' + (£- V Po ) = p' + (€ • e r )^ =p' + (10) 

where e r is the radial unit vector. 

In terms of the Eulerian perturbations and the displacement vector, £ the 
linearized mass, momentum and energy equations can be expressed in the 
following form: 

p' + V( P oO=0, (11) 

dv 

Po-t^ = -VP' - g e r p' + p V<Z>', (12) 

P +Cr^ =C Q (p +^ r — ), (13) 

V 2 <2>' = 4nGp'. (14) 

The equation of solar oscillations can be further simplified by neglecting 
the perturbations of the gravitational potential, which give relatively small 
corrections to theoretical oscillation frequencies. This is so-called Cowling 
approximation: <P' = 0. 

Now, we consider the linearized equations in the spherical coordinate sys- 
tem, r, 9, </>. In this system, the displacement vector has the following form: 

i = £ r e r + S,eee + = £ r e r + (15) 

where £h = Ce e 6» + £c/>e<p is the horizontal component of displacement. Also, 
we use the equation for divergence of the displacement (called dilatation): 

V* ee d^ =\^~(r%) + -^^(sin^) + -^-f-t = 
r z or r sin o0 r sin oq> 

= ^%) + \v h i h . (16) 
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We consider periodic perturbations with frequency ui: £ oc exp(zwi), .... Here, 
uj is the angular frequency measured in rad/sec; it relates to the cyclic 
frequency, v, which measures the number of oscillation cycles per sec, as: 
UJ = 2ttv. 

Then, in the Cowling approximation, we obtained the following system of 
the linearized equations (omitting subscript for unperturbed variables): 

QP> 

-uj 2 p^ = -—+gp', (18) 

-uj 2 pti h = --V h P', (19) 
r 

p> = \p> + ^ r , (20) 
c 2 g 



where 

N 2 = n( 

K -fP dr P dr / 

is the Briint-Vdisdld (or buoyancy) frequency. 

For the boundary conditions, we assume that the solution is regular at 
the Sun's center. This correspond to the zero displacement, £ r = at r = 0, 
for all oscillation modes except of the dipole modes of angular degree 1 = 1. 
In the dipole-mode oscillations the center of a star oscillates (but not the 
center of mass), and the boundary condition at the center is replaced by a 
regularity condition. At the surface, we assume that the Lagrangian pressure 
perturbation is zero: SP = at r = R. This is equivalent to the absence 
of external forces. Also, we assume that the solution is regular at the poles 

9 = 0,77. 

We seek a solution of Eqs (17-20) by separation of the radial and angular 
variables in the following form: 

p'(r,9^)=pi(r)-f(9^), (22) 

P'(r,9^)=P'(r)-f(9^), (23) 

Zr(r,6,<t>)=Zr(r)-M4>), (24) 

€h(r,M)=&(r)V h /(M). (25) 

Then, in the continuity equation: 
1 d 



f(6,cj>) + ^ h V 2 h f = 0. (26) 



the radial and angular variables can be separated if 



VLf = af, 



(27) 



26 Alexander G. Kosovichev 

where a is a constant. 

It is well-known that this equation has a non-zero solution regular at the 
poles (9 — 0, 7r) only when 

a = -1(1 + 1), (28) 
where I is an integer. This non-zero solution is: 

f(9,ct>)=Y l m (e,4>)<xPr(6)e im 'l>, (29) 

where PJ n (9) is the associated Legendre function of angular degree I and 
order to. 

Then, the continuity equation for the radial dependence of the Eulerian 
density perturbation, p'(r), takes the form: 

p' + ^|:(^)-^^ = o. (30) 

The horizontal component of displacement £/, can be determined from the 
horizontal component of the momentum equation: 

-Lo 2 pUr) = -- r P\r), (31) 

or 

a = -4-p*- (32) 

Substituting this into the continuity equation (30) we get: 

di r dp 2 ^ P' pN 2 ^ L 2 . n 

p-r + ih-f + -pir + — + —ir - -^r p = o, 33 

dr dr r c 2 g r 2 uj 2 p 

where we define L 2 = 1(1 + 1). 

Using the hydrostatic equation for the background (unperturbed) state, 
dP/dr = —gp, we finally obtain: 



or 



where 



L 2 c 2 ' 


)4 


r 2 u) 2 1 


/ pc 2 


sf\ 


P' 


oj 2 J 


pc 2 



T 2 r 2 

Sf = (36) 



is the Lamb frequency. 

Similarly, for the momentum equation we obtain: 

dP ' , 9 D / , /Ar2 ,2 



P' + (N - LU 2 )pCr = 0. (37) 
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The inner boundary condition at the Sun's center is: 

Cr = 0, (38) 

or a regularity condition for 1 = 1. 

The outer boundary condition at the surface (r = R) is: 

5P = P' + ^f-Cr = 0. (39) 
ar 

Applying the hydrostatic equation, we get: 

P' - gp(; r = 0. (40) 

Using the horizontal component of the momentum equation: P' = oj pr£h, 
the outer boundary condition (40) can be written in the following form: 

^ = 4-, (41) 

that is the ratio of the horizontal and radial components of displacement is 
inverse proportional to the squared oscillation frequency. However, observa- 
tions show that this relation is only approximate, presumably, because of the 
external force caused by the solar atmosphere. 



p-modc (7=20, m=16, «=14) g-mode (1=5, »i=3, «=6) 




Fig. 15. Eigenfunctions (42) of two normal oscillation modes of the Sun: a) p-mode 
of angular degree I = 20, angular degree m = 16, and radial order n = 16, b) g- 
mode of I = 5, m = 3, and n = 5. Red and blue-green colors correspond to positive 
and negative values. 



Equations (35) and (37) with boundary conditions (38)- (40) constitute an 
eigenvalue problem for solar oscillation modes. This eigenvalue problem can 
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be solved numerically for any solar or stellar model. The solution gives the 
frequencies, cj n i, and the radial eigenfunctions, £r™^(r) and P'^ n ' l \r), of the 
normal modes. 

The radial eigenfunctions multiplied by the angular eigenfunctions (22)- 
(25) represented by the spherical harmonics (29) give three-dimensional os- 
cillation eigenfunctions of the normal modes, e.g.: 

Ur,e,M=e 7 n ' l \r)Yr{e,<t>). (42) 

Examples of such two eigenfunctions for p- and g-modes are shown in Fig. 15. 
It illustrates the typical behavior of the modes: the p-modes are concentrated 
(have the strongest amplitude) in the outer layers of the Sun, and g-modes 
are mostly confined in the central region. 

4.2 JWKB solution 

The basic properties of the oscillation modes can be investigated analyti- 
cally using an asymptotic approximation. In this approximation, we assume 
that only density p(r) varies significantly among the solar properties in the 
oscillation equations, and seek for an oscillatory solution in the JWKB form: 

Cr = Ap- 1/2 e lkrr , (43) 
P' = Bp l ' 2 e lkr \ (44) 

where the radial wavenumber k r is a slowly varying function of r; A and B 
are constants. 

Then, substituting these in Eqs (35) and (37) we obtain: 



(45) 



"'" = -Bp" 2 (-it,--!-)e"- r , (46) 



dr ' V H 

where 

dr 

is the density scale height. 

From (45-46) we get a linear system for the constant, A, and B: 

-^ + £^-^+1(1-5)5 = 0, (48) 

-ik r -^-)B+^B + (N 2 -Lj 2 )A = 0. (49) 
HI c 
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It has a non-zero solution when the determinant is equal zero, that is when 



where 



2H 



(50) 



(51) 



is the acoustic cut-off frequency. Here, we used the relation: iV 2 = g/H 
9 2 /c\ 
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Fig. 16. Buoyancy (Brunt- Vaisala) frequency N (thick curve), acoustic cut-off fre- 
quency, lu c (thin curve) and Lamb frequency Si for 1=1, 5, 20, 50, and 100 (dashed 
curves) vs. fractional radius r/R for a standard solar model. The horizontal lines 
with arrows indicate the trapping regions for a g mode with frequency v = 0.2 
mHz, and for a sample of five p modes: I = 1, v = 1 mHz; I = 5, v = 2 mHz; I — 20, 
v = 3 mHz; I = 50, v = 4 mHz; / = 100, v = 5 mHz. 



The frequencies of solar modes depend on the sound speed, c, and three 
characteristic frequencies: acoustic cut-off frequency, oj c (51), Lamb frequency, 
Si (36), and Brunt- Vaisala frequency, N (21). These frequencies calculated 
for a standard solar model are shown in Fig. 16. The acoustic cut-off and 
Briint- Vaisala frequencies depend only on the solar structure, but the Lamb 
frequency depends also on the mode angular degree, I. This diagram is very 
useful for determining the regions of mode propagation. The waves propagate 
in the regions where the radial wavenumber is real, that fc 2 > 0. If fc 2 < 
then the waves exponentially decay with distance (become 'evanescent'). The 
characteristic frequencies define the boundaries of the propagation regions, 
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also called the wave turning points. The region of propagation for p- and 
g-modes are indicated in Fig. 16, and are discussed in the following sections. 
We define a horizontal wavenumber as 

k h = -, (52) 



where L = y/l(l + 1). This definition follows from the angular part of the 
wave equation (27): 

^irr + ll ^Yr = o, (53) 

where Vh is the horizontal component of gradient. It can be rewritten in terms 
of a horizontal wavenumber, k h , ^V/^" 1 + fc 2 !?" = if k 2 = 1(1 + l)/r 2 . 
In term of kh the Lamb frequency is Si — k h c, and Eq. 50 takes the form: 

«(£-')• (54 » 

The frequencies of normal modes are determined for the Borh quantization 
rule (resonant condition): 

r 2 

k r dr = ir(n + a), (55) 

i 

where n and T2 are the radii of the inner and outer turning points where 
fc r =0, n is a radial order -integer number, and a is a phase shift which depends 
on properties of the reflecting boundaries. 



4.3 Dispersion relations for p- and g-modes 

For high-frequency oscillations, when lu 2 >> N 2 , the dispersion relation (50)- 
(54) can be written as: 

,,2 , ,2 c2 , ,2 , ,2 
fc r - c2 - c2 - c2 - (bOj 

Then, we obtain: 

u? = u? c + (k 2 + k 2 h )c 2 = uj 2 c + k 2 c 2 . (57) 

This is a dispersion relation for acoustic (p) modes, uo c is the acoustic cut-off 
frequency. The wave with frequencies less than uj c (or wavelength A > AnH) 
do not propagate. These waves exponentially decay, and called 'evanescent'. 
For low- frequency perturbations, when uj 2 « Sf, one gets: 

fc r 2 = #^ 2 -^) = %iV 2 - W 2 ), (58) 
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h 2 N 2 

u 2 = = N 2 cos 2 9, (59) 

where 9 is the angle between the wavevector, k, and horizontal surface. 

These waves are called internal gravity waves or g-modes. They propagate 
mostly horizontally, and only if lu 2 < N 2 . The frequency of the internal 
gravity waves does not depend on the wavenumbcr, but on the direction of 
propagation. These waves are evanescent if lo 2 > N 2 . 



4.4 Frequencies of p- and g-modes 

Now, we use the Borh quantization rule (55) and the dispersion relations for 
the p- and g-modes (57-58) to derive the mode frequencies. 

p-modes: The modes propagate in the region where k 2 > 0; and the radii of 
the turning points, r\ and r 2 , are determined from the relation k 2 = 0: 



2 _ , ,2 



< + -o- = °- (6°) 



The acoustic cut-off is only significant near the Sun's surface. The lower 
turning point is located in the interior where io c « uj (Fig. 16. Then, at the 
lower turning point, r = r\\ oj w Lc/r, or 

(61) 

represents the equation for the radius of the lower turning point, r\. The 
upper turning point is determined by the acoustic frequency term: uj c (r 2 ) ~ u. 
Since w c (r) is a steep function of r near the surface, then 

r 2 w R. (62) 

The p-mode propagation region is illustrated in Fig. 16 Thus, the resonant 
condition for the p-modes is: 



L 2 

^ T dr = ir(n + a) (63) 



In the case, of the low-degree "global" modes, for which I « n, the lower 
turning point is almost at the center, n w 0, and we obtain [17]: 
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This relation shows is the spectrum of low-degree p-modes is approximately 
equidistant with the frequency spacing: 



Av = 4 / — . (65) 




This corresponds very well to the observational power spectrum shown in 
Fig. 4. According to this relation, the frequencies of mode pairs, (n, I) and 
(n — 1,1 + 2), coincide. However, calculations to the second-order shows that 
the frequencies in these pairs are separated by the amount [98,99]: 

Av f R dc dr 

Sv nl = v nl - v n _ hl+2 « -{Al + 6)^- ^ --. (66) 

This is so-called "small separation". For the Sun, Av rts 136^iHz, and 5v rts 
9^iHz. The I — v for the p-modes is illustrated in Fig. 17. 

g-modes: The turning points, k r — 0, are determined from equation (58): 

N(r) = lo. (67) 

In the propagation region, k r > 0, (see Fig. 16), far from the turning points 

(A^ >>u;): 

k r w — . (68) 

ruj 
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Then, from the resonant condition: 



N—=ir{n + a). (69) 



we find an asymptotic formula for the g-mode frequencies: 

WRj Jn r _ 

7r(n + a) 



(70) 



It follows that for a given I value the oscillation periods form a regular equally 
spaced pattern: 

2n n(n + a) 

~ = L r N— ' 

Jri r 

The distribution of numerically calculated g-mode periods is shown in Fig. 4.4. 



p= — = "y^^> (7i) 




i 2 3 4 5 G 7 8 9 10 
angular degree, i 



Fig. 18. Periods of solar oscillation modes in the angular degree range, I = — 10. 
Labels <7i-<?6 mark g-modes of the radial order n = 1, . . . , 6 



4.5 Asymptotic ray-path approximation 

The asymptotic approximation provides an important representation of solar 
oscillations in terms of the ray theory. Consider the wave path equation in 
the ray approximation: 

* = (72) 

dt dk y ' 
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Then, the radial and angular components of this equation are: 

- - — (73) 

r dt = w h (74) 

Using the dispersion relation for acoustic (p) modes: 

o, 2 = c 2 (fc 2 + fc 2 ), (75) 

in which we neglected the ui c term. (It can be neglected everywhere except 
near the upper turning point, R), we get 

(IT 

dt= m . (76) 

Thus, is the travel time from the lower turning point to the surface. 

The equation for the acoustic ray path is given by the ratio of equations 
(74) and (73): 

dB = ( duj \ i ( duj_\ = kh (77) 
T dr \dk h )'\dk r ) kr' 

or 

dB = k ± = L/r 

dr k r y/uj 2 /c 2 - L 2 /r 2 

For any given values of co and /, and initial coordinates, r and 8, this equation 
gives trajectories of ray paths of p-modes inside the Sun. The ray paths 
calculated for two solar p-modes are shown in Fig. 19a. They illustrate an 
important property that the acoustic waves excited by a source near the solar 
surface travel into the interior and come back to surface. The distance, A, 
between the surface points for one skip can be calculated as the integral: 



A = 2f R dO = 2[ R L ' r ir = 2 f ± ir 

J ri J ri ^ 2 /c 2 - L 2 /r 2 J ri V^ 2 / L2 - c2 / r2 

(79) 

The corresponding travel time is calculated by integrating equation (73): 

* * (80) 



= 2 r dt =r — 

J ri Jrx c(l-fc 2 C 2 /w 2 ) 1/2 Jr 



1/2- 



(l-k 2 h c 2 /uj 2 y /z J ri c(l- L 2 c 2 /r 2 L0 2 ) 

These equations give a time- distance relation, t~A, for acoustic waves travel- 
ing between two surface points through the solar interior. The ray representa- 
tion of the solar modes and the time-distance relation provided a motivation 
for developing time- distance helio seismology (Sec. 7), a local helioseismology 
method [48]. 
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Fig. 19. Ray paths for a) two solar p-modes of angular degree I = 2, frequency 
v — 1429.4 ^iHz (thick curve), and I = 100, v = 3357.5 piHz (thin curve); b) g-mode 
of I = 5, v = 192.6 /iHz (the dotted curve indicates the base of the convection 
zone). The lower turning points, n of the p-modes are shown by arrows. The upper 
turning points of these modes are close to the surface and not shown. For the g- 
mode, the upper turning point, r2, is shown by arrow. The inner turning point is 
close to the center and not shown. 



The ray paths for g-modes are calculated similarly. For the g-modes, the 
dispersion relation is: 

k 2 N 2 

^ = • ( 81 ) 

k r + k h 

Then, the corresponding ray path equation: 



d0 K N 2 , 

r Tr=-i=-y^- 1 - (82) 

The solution for a g-mode of I = 5, v = 192.6/iHz is shown in Fig. 19b. 
Note that the g-mode travels mostly in the central region. Therefore, the 
frequencies of g-modes are mostly sensitive to the central conditions. 



4.6 Duvall's law 

The solar p-modes, observed in the period range of 3-8 minutes, can be 
considered as high-frequency modes and described by the asymptotic theory 
quite accurately. Consider the resonant condition (63) for p-modes: 



f R fuj 2 i 2 \ 1/2 , 



(83) 
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Dividing both sides by lj we get: 



ir(n + a) 



LO 



(84) 



Since the lower integral limit, n depends only on the ratio L/co, then the 
whole left-hand side is a function of only one parameter, L/co, that is: 



F 



L\ 7r(n + a) 



(85) 



This relation represents so-called Duvall's law [100]. It means that a 2D 
dispersion relation lo = co(n, I) is reduced to the ID relation between two 
ratios L/lo and (n+a)/u>. With an appropriate choice of parameter a (e.g. 1.5) 
these ratios can easily calculated from a table of observed solar frequencies. 
An example of such calculations shown in Fig. 20) illustrates that the Duvall's 
law holds quite well for the observed solar modes. The short bottom branch 
that separates from the main curve correspond to f-modes. 
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Fig. 20. The observed Duvall's law relation for modes of I = — 250. 



4.7 Asymptotic sound-speed inversion 

The Duvall's law demonstrate that the asymptotic theory provides a rather 
accurate description of the observed solar p-modes. Thus, it can be used for 
solving the inverse problem of helioseismology - determination of the internal 
properties from the observed frequencies. Theoretically, the internal structure 
of the Sun is described by the stellar evolution theory [101]. This theory cal- 
culates the thermodynamic structure of the Sun during the evolution on the 
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Main Sequence. The evolutionary model of the current age, w 4.6 x 10 9 years, 
is called the standard solar model. Helioseismology provides estimates of the 
interior properties, such as the sound-speed profiles, that can be compared 
with the predictions of the standard model. 

Our goal is to find corrections to a solar model from the observed fre- 
quency differences between the Sun and the model using the asymptotic 
formula for the Duvall's law [102]. 

We consider a small perturbation of the sound-speed, c — > c + Ac, and the 
corresponding perturbation of frequency: ui — > to + Auj. Then, from equation 
(84) we obtain: 



R r (u + Auj) 2 L 2 
(c + Ac) 2 r 2 



1/2 

dr = ir(n + a). (86) 



Expanding this in terms of Ac/c and Aui/ui and keeping only the first-order 
terms we get: 



Aw f R dr f 

~ Jr t c(l- L 2 c 2 /r 2 tu 2 ) 1/2 ~ Jr 



R Ac dr 



c(l- L 2 c 2 /r 2 uj 2 ) 1/2 Jr t c c{\- L 2 c 2 /r 2 uo 2 ) 1/2 ' 
If we introduce a new variable: 



(87) 



-r 

J r t 



In c(l-L 2 c 2 /r 2 Lj 2 ) 1/2 

then 



(88) 



Auj 1 f R Ac dr 



w T Jrt c c(l- L 2 c 2 /r 2 io 2 ) 1/2 ' 

This equation has a simple physical interpretation: T is the travel time of 
acoustic waves to travel along the acoustic ray path between the lower and 
upper turning points (Fig. 19). The right-hand side integral is an average of 
the sound-speed perturbations along this ray path (compare with Eq.(80)). 

Equation (89) can be reduced to the Abel integral equation by making a 
substitution of variables. The new variables are: 

x=^, (90) 

V = ^5> ( 91 ) 

where £ is a measured quantity, and y is associated with the sound-speed 
distribution of an unperturbed solar model. 
Then, we obtain an equation for x and y: 



(92) 
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where 

f(y) = 



C 2j/ 3/ 2 f^g£ + l 

\ d log r 



To solve for f(y) we multiply both sides of Eq.(12) by dx/^/z — x and 
integrate with respect to x from to z: 

F{x)dx _ f z dx [ x f(y)dy _ 



o y/z - x Jo yjz - x Jo \Jx-y 

dx 



r f(y)dy r 

JO Jy 



\/{z- x)(x - y) 



Here we changed the order of integration. 
Note that 

dx 

>y \f{z-x){x-y) 

then 



Jy 



Jo y/Z-X Jo 

Differentiating with respect to x, we obtain the final solution: 

m-iff ™*. (93) 

n dx Jo \JZ — X 

Then, from f(y) we find the sound-speed correction Ac/c. 

This method based on linearization of the asymptotic Abel integral is 
called "differential asymptotic sound-speed inversion" [102]. It provides esti- 
mates of the sound-speed deviations from a reference solar model. 

Alternatively, the sound-speed profile inside the Sun can be found from 
a implicit solution of the Abel obtained by differentiating the Duvall's law 
equation (84) with respect to variable y — L/uj. Then, this equation can be 
solved analytically. The solution provides an implicit relationship between 
the solar radius and sound speed [103]: 

rR/c, jp / 2\ -1/2 

Hr/R)= y 2- ) d y, (94) 



r/c dy 

where c s is the sound speed at the solar surface r — R. The calculation 
of the derivative, dF/dy, is essentially differentiation of a smooth function 
approximating the Duvall's law, that is differentiating 7r(n+a)/w with respect 
to L/co. Both of these quantities are obtained from the observed frequency 
table, co(n, I). 
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The first inversion results using this approach was published by Christensen- 
Dalsgaard et al [102]. These technique can be generalized by including the 
Brunt- Vaisala frequency term in the p-mode dispersion relation, and also 
taking into account the frequency dependence of the phase shift, a [36]. The 
results show that this inversion procedure provides a good agreement with the 
solar models, used for testing, except the central core, where the asymptotic 
and Cowling approximations become inaccurate. 

3 I . . ■ 1 0.04 1 1 1 1 1 1 1 1 1 1 1 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

a) r/R b) r/R 

Fig. 21. a) Result of the asymptotic sound inversion (solid curve) [104] for the 
p-mode frequencies [105]. It confirmed the standard solar model (model 1) [106] 
(dots). The large discrepancy in the central region is due inaccuracy of the data 
and the asymptotic approximation, b) The relative difference in the squared sound 
speed between the asymptotic inversions of the observed and theoretical frequencies. 



Figure 21 shows the inversion results [104] for the p-mode frequencies 
measured by Duvall ct al. [105]. The deviation of the sound speed from a 
standard solar model is about 1%. Later, the agreement between the solar 
model and and the helioseismic inversions was improved by using more precise 
opacity tables and including element diffusion in the model calculations [101]. 
Also, a more accurate inversion method was developed by using a perturba- 
tion theory based on a variational principle for the normal mode frequencies 
(Sec. 5). 

4.8 Surface gravity waves (f-mode) 

The surface gravity (f-mode) waves are similar in nature to the surface ocean 
waves. They are driven by the buoyancy force, and exist because of the sharp 
density decrease at the solar surface. These waves are missing in the JWKB 
solution. These waves propagate at the surface boundary where Lagrangian 
pressure perturbation SP <~ 0. 

To investigate these waves we consider the oscillation equations in terms of 
5P by making use of the relation between Eulerian and Lagrangian variables 
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(10): 

P' = 5P + gp£ r . 

The oscillation equations (35) and (37) in terms of £ r and SP are: 

& + 1 _ 73 = °> ( 95 ) 



dr w 2 r 2 \ to 2 r 2 J pc 2 

dSP + ^5P-^ r = 0, (96) 



dr uj 2 r 2 r 
where 

i 2 .9 



/« 5^. (97) 



These equations have a peculiar solution: 

5P = 0, / = 0. 

For this solution: 



u? = ^ = k h g (98) 



-dispersion relation for f-mode. 
The eigenfunction equation: 



^ - "€r = (99) 

has a solution 

£ r oc e fc * (»-*) (100) 

exponentially decaying with depth. 

These waves are similar in nature to water waves which have the same 
dispersion relation: uj = gkh- The f-mode waves are incompressible: Vv = 0. 
These waves are not sensitive to the sound speed but are sensitive to the 
density gradient at the solar surface. They are used for measurements of the 
'seismic radius' of the Sun. 



4.9 The seismic radius 

The frequencies of f-modes: 

2 , GM L GM 
U =^ S ^"5 SL ^- (101) 

If the frequencies are determined in observations for given I, then we can 
define the 'seismic radius', R, as 
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The procedure of measuring the solar seismic radius is simple [107]. The 
lower curve in Figure 22a shows the relative difference between the f-modc 
frequencies of I = 88—250 calculated for a standard solar model (Model S) and 
the frequencies obtained from the SOHO/MDI observations. This difference 
shows that the model frequencies are systematically, by ~ 6.6 x 10 -4 , lower 
than the observed frequencies. Then from equation (101): 



AR 
R 



lAv 



4.4 x 10" 



(103) 



This means that the seismic radius is approximately equal to 695.68 Mm, 
which is about 0.3 Mm less than the standard radius, 695.99 Mm, used for 
calibrating the model calculation. This radius is usually measured astromet- 
rically as a position of the inflection point in the solar limb profile. However, 
in the model calculations it is considered as a height where the optical depth 
of continuum radiation is equal 1. The difference between this height and the 
height of the inflection point can explain the discrepancy between the model 
and seismic radius. 




-lxltr 



695.0 695.5 
radius (Mm) 



Fig. 22. a) Relative differences between the f-mode frequencies of I = 88 — 250 
computed for a standard solar model (Model S) and the observed frequencies. 
The 'seismic model' frequencies are obtained by scaling the frequencies of model 
S with factor 1.00066 which corresponds to scaling down the model radius with 
(1.00066) 2/3 » 1.00044. The error bars are 3<r error estimates of the observed fre- 
quencies, b) Density as a function of radius near the surface for the standard and 
seismic models. The star indicates the photospheric radius. The diamond shows the 
seismic radius, 695.68 Mm. 



Figure 22b illustrates the density profiles in the standard solar model 
(model S [101]) and a 'seismic' model, calibrated to the seismic radius. The 
f-mode frequencies of the seismic model match the observations. 
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1.1 1.2 1.3 1.4 1.5 1.6 1.7 

v (iilHz) 

Fig. 23. Average relative frequency differences in f-mode (8v/u) as a function of 
{u), average frequencies binned every 20 fiHz. The reference year is 1996. 



Since the f-mode frequencies provide an accurate estimate of the seismic 
radius, then it is interesting to investigate the variations of the solar radius 
during the solar activity cycle, which are quite important for understanding 
physical mechanisms of solar variability (e.g. [108]). Figure 23 shows the f- 
mode frequency variations during the solar cycle 24, in 1997-2004, relative to 
the f-mode frequencies observed in 1996 during the solar minimum [109]. 

The results show a systematic increase of the f-mode frequency with the 
increased solar activity, which means a decrease of the seismic radius. How- 
ever, the variations of the f-mode frequencies are not constant as this is 
expected from equation (103 for a simple homologous change of the solar 
structure. A detailed investigation of these variations showed that the fre- 
quency dependence can be explained if the variations of the solar structure 
are not homologous and if the deeper subsurface layers expand but the shal- 
lower layers shrink with the increased solar activity [109,110]. 

5 General helioseismic inverse problem 

In the asymptotic (high-frequency of short-wavelength) approximation (84), 
the oscillation frequencies depend only on the sound-speed profile. This de- 
pendence is expressed in terms of the Abel integral equation (89), which can 
be solved analytically. 

In the general case, the relation between the frequencies and internal prop- 
erties is more complicated, the frequencies depend not only on the sound 
speed, but also on other internal properties, and there is no analytical so- 
lution. Generally, the frequencies determined from the oscillation equations 
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(35) and (37) depend on the density, p(r), the pressure, P(r), and the adia- 
batic exponent, 7(r). However, p and P are not independent, and related to 
each other through the hydrostatic equation: 

^ = -™ < 104 > 

where g = Gm/r 2 , m = An J Q r pr' 2 dr' . Therefore, only two thermodynamic 
(hydrostatic) properties of the Sun are independent, e.g. pairs of (p, 7), (P, 7), 
or their combinations: (P/p, 7), (c 2 ,7), (c 2 ,p) etc. 

The general inverse problem of helioseismology is formulated in terms of 
small corrections to the standard solar model because the differences between 
the Sun and the standard model are typically 1% or less. When necessary 
the corrections can be applied repeatedly, using an iterative procedure. 



5.1 Variational principle 

We consider the oscillation equations as a formal operator equation in terms 
of the vector displacement, £: 

u 2 t = £(£), (105) 

where C in the general case is an integro-differential operator. If we multiply 
this equation by and integrate over the mass of the Sun we get: 

u 2 I pC ■ £dV = / £* ■ CipdV, (106) 
Jv Jv 

where p is the model density, V is the solar volume. 

Then, the oscillation frequencies can be determined as a ratio of two 
integrals: 

co 2 = r -z ; JTA ■ 107 

The frequencies are expressed in terms of eigenfunctions ^ and the solar 
properties properties represented by coefficients of the operator C. For small 
perturbations of solar parameters the frequency change will depend on these 
perturbations and the corresponding perturbations of the eigenfunctions, e.g. 

5lu 2 = &[5p,5f,5€\. (108) 

The variational principle states that the perturbation of the eigenfunc- 
tions constitute second-order corrections, that is to the first-order approx- 
imation the frequency variations depend only on variations of the model 
properties: 

Sco 2 w V[8p,5i\. (109) 
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The variational principle allows us to neglect the perturbation of the eigen- 
functions in the first-order perturbation theory. This was first established by 
Rayleigh. Thus, equation (107) is called the Rayleigh's Quotient, and the 
variational principle is called the Rayleigh's Principle. The original formula- 
tion of this principle is: for an oscillatory system the averaged over period 
kinetic energy is equal the averaged potential energy. In our case, the left- 
hand side of equation (106) is proportional to the mean kinetic energy, and 
the right-hand side is proportional to the potential energy of solar oscillations. 



5.2 Perturbation theory 

We consider a small perturbation of the operator C caused by variations of 
the solar structure properties: 

£(€) = £o(0 + £i(€)- 

Then, the corresponding frequency perturbations are determined from the 
following equation: 



Slj 



or 

5uj 1 



: __j v e-c^ P dv 



where 



! v pi*-idv ' 

r • c^pdv, (no) 



v 



I = 



[ p£*-£dV (ill) 
Jv 



is so-called mode inertia or mode mass. The mode energy is E = Iui^a , where 
a is the amplitude of the surface displacement. The mode eigenfunctions are 
usually normalized such that £ r (-R) = 1. 

Using explicit formulations for operator L\ Eq. 110 can be reduced to a 
system of integral equations for a chosen pair of independent variables [111— 
114], e.g. for (p,7) 



[ R K (n,D S _P dr+ j R M)h d ^ (n2) 

Jo P Jo 7 



where Kp U j\r) and K^J\r) are sensitivity (or 'seismic') kernels. They are 
calculated using the initial solar model parameters, p , P , 7, and the oscil- 
lation eigenfunctions for these model, £. 



5.3 Kernel transformations 

The sensitivity kernels for various pairs of solar parameters can be obtained 
by using the relations among these parameters, which follows from the equa- 
tions of solar structure ('stellar evolution theory'). 
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A general procedure for calculating the sensitivity kernels developed by 
Kosovichev [114] can be illustrated in an operator form. Consider two pairs 
of solar variables, X and Y, e.g. 

Sp SY\ _ (6u 5Y_ 
p ' 7 J ' V u ' Y 

where u = P/p, Y is the helium abundance. 

The linearized structure equations (the hydrostatic equilibrium equation 
and the equation of state) that relates these variables can be written sym- 
bolically: 

AX = Y. (113) 

Let Kx and Ky be the sensitivity kernels for X and Y, then the fre- 
quency perturbation is: 

' R K x -Xdr=(K x -X), (114) 



w Jo 

where < • > denotes the inner product. Similarly, 

— = (K Y -Y). (115) 

Then from equations (114) and (115) we obtain the following relation: 

(K Y ■ Y) = (K Y ■ AX) = (A*K Y ■ X) , (116) 

where A* is an adjoint operator. This operator is adjoint to the stellar struc- 
ture operator, A. The second part of equation (116) represent a formal defi- 
nition of this operator. 

From Eq.(114) and (116) we get: 

(A*K Y ■ X) = (K x ■ X) . 

This equation is valid for any X only if 

A*K Y = K x . (117) 

That means that the equation for the sensitivity kernels is adjoint to the 
stellar structure equations. The explicit formulation of the adjoint equations 
for the sensitivity kernels for various pairs of variables is given in [114]. 

Examples of the sensitivity kernels for solar properties are shown in Fig- 
ures 24. Figure 25 illustrates the difference in sensitivities of the p- and g- 
modes. The frequencies of solar p-modes are mostly sensitive to properties of 
the outer layers of the Sun while the frequencies of g-modes have the greatest 
sensitivity to the parameters of the solar core. 
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Fig. 24. Sensitivity kernels for the acoustic mode of the angular degree, I— 10, 
and the radial order, n—6. K Ptl is the kernel for density, p, at constant adiabatic 
exponent, 7; K c i p is the kernel for the squared sound speed, c 2 , at constant p; K u> y 
is the kernel for function u, - the ratio pressure, p, to density at constant helium 
abundance, Y; and K a* , 7 is the kernel for the parameter of convective stability, 
A* — rN 2 /g, at constant 7. 



5.4 Solution of inverse problem 

The variation formulation provides us with a system integral equations (112) 
for a set of observed mode frequencies. Typically, the number of observed fre- 
quencies, N ~ 2000. Thus, we have a problem of determining two functions 
from this finite set of measurements. In general, it is impossible to determine 
these functions precisely. We can always find some rapidly oscillating func- 
tions, f(r), that being added to the unknowns, 5p/p and #7/7, do not change 
the values of the integrals, e.g. 

/\M(r)/(r)dr = 0. 
Jo 

Such problems without a unique solution are called " ill-posed" . The gen- 
eral approach is to find a smooth solution that satisfies the integral equations 
(112) by applying some smoothness constraints to the unknown functions. 
This is called a regularization procedure. 

There are two basic methods for the helioseismic inverse problem: 
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Fig. 25. Sensitivity kernels for p- and g-modes for u — P/ p and helium abundance 
Y. 



1. Optimally Localized Averages (OLA) method - (Backus-Gilbert method) 
[115] 

2. Regularized Least-Squares (RLS) method - (Tikhonov method) [116] 
5.5 Optimally localized averages method 

The idea of the OLA method is to find a linear combination of data such 
as the corresponding linear combination of the sensitivity kernels for one 
unknown has an isolated peak at a given radial point, ro, (resembling a 5- 
function), and the combination for the other unknown is close to zero. Then, 
this linear combination provides an estimate for the first unknown at ro. 

Indeed, consider a linear combination of (112) with some unknown coef- 
ficient a^™'^: 

y a (n,i) s _^l == f R y a MM)Sp dr+ [ R y a (n,D K (ni)h dr ^ 

w (n,I) J P<1 p J ^ 7 ' P 7 

(118) 

If in the first term the linear combination of the kernels is close to a 5-function 
at r = r , that is 

^fl<"%W>(r)^( r -r„), (119) 
and the linear combination in the second term vanishes: 



^a("%M(r)~0, 



(120) 
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then equation (118) gives an estimate of the density perturbation, Sp/p, at 
r = r : 




(121) 



Of course, the coefficients, a( n - l \ of equation (121) must be calculated from 
conditions (119) and (120) for various target radii r . 
The functions, 

^a< n MKjtf\r) = A(r ,r), (122) 

^a("%W»(r)E5(ro,r), (123) 

are called the averaging kernels. They play a fundamental role in the he- 
lioseismic inverse theory for determining the resolving power of helioseismic 
data. 

The coefficients, a n ' 1 , are determined my minimizing a quadratic form : 

f R 

M(r ,A,a,(3)= J{r ,r) [A(r ,r)f dr + (124) 
Jo 

+(3 [ [B(r ,r)} 2 dr + a y2E n , l;n , M a n < l a n '' 1 ' , 

where function J(r , r) = 12(r — r ) 2 provides a localization of the averaging 
kernels A(r,r ) at r = r , E n j. n >^ is a covariance matrix of observational 
errors, a and ft are regularization parameters. The first integral in cq. (125) 
represents the Backus-Gilbert criterion of J-ness for A(r ,r); the second term 
minimizes the contribution from B(ro,r), thus, effectively eliminating the 
second unknown function, (#7/7 in this case); and the last term minimizes the 
errors. A practical minimization algorithm is presented in [114]. An example 
of the averaging kernels is shown in Fig. 26 



5.6 Inversion results for solar structure 

As an example, consider the results of inversion of the recent data obtained 
from the MDI instrument on board the SOHO space observatory. The data 
represent 2176 frequencies of solar oscillations of the angular degree, I, from 
to 250. These frequencies were obtained by fitting peaks in the oscillation 
power spectra from a 360-day observing run, between May 1, 1996 and April 
25, 1997. 

Figure 27 shows the relative frequency difference, Suj/uj, between the 
observed frequencies and the corresponding frequencies calculated for the 
standard model S [101]. The frequency difference is scaled with a factor 
Q = I(co)/Io(co), where I{uj) is the mode inertia, and Iq{oj) is the mode 
inertia of radial modes (/ = 0), calculated at the same frequency. 
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Fig. 26. A sample of the optimally localized averaging kernels for the structure 
function, u, the ratio of pressure, P, to density, p,u — P/p. The second, eliminated, 
parameter in these kernels is the helium abundance, Y. 
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Fig. 27. The relative frequency difference, scaled with the relative mode inertia 
factor, Q — I/I (111), between the Sun and the standard solar model. 



This scaled frequency difference depends mainly on the frequency alone 
meaning that most of the difference between the Sun and the reference solar 
model is in the near-surface layers. Physically, this follows from the fact that 
the p-modes of different I behave similarly near the surface where they prop- 
agate almost vertically. This behavior is illustrated by the p-mode ray paths 
in Fig. 19a, which become almost radial near the surface. In the inversion 
procedure, this frequency dependence is eliminated by adding an additional 
"surface term" in equation (112) [114]. However, there is also a significant 
scatter along the general frequency trend. This scatter is due to the variations 
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of the structure in the deep interior, and it is the basic task of the inversion 
methods to uncover the variations. 
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Fig. 28. The results of test inversions (points with the error bars, connected with 
dashed curves) of frequency differences between two solar models for the squared 
sound speed, c 2 , the adiabatic exponent, 7, the density, p, and the parameter of 
convective stability, A*. The solid curves show the actual differences between the 
two models. Random Gaussian noise was added to the frequencies of a test solar 
model. The vertical bars show the formal error estimates, the horizontal bars show 
the characteristic width of the localized averaging kernels. The central points of the 
averages are plotted at the centers of gravity of the averaging kernels. 



First, we test the inversion procedure by considering the frequency differ- 
ence for two solar models and trying to recover the differences between model 
properties. Results of the test inversion (Fig. 28) show good agreement with 
the actual differences. However, the sharp variations, like a peak in the pa- 
rameter of convective stability, A* = rN 2 /g, at the base of the convection 
zone, are smoothed. Also, the inner 5% of the Sun and the subsurface layers 
(outer 2-3%) are not resolved. 

Then, we apply this procedure to the real solar data. The results (Fig. 29) 
show that the differences between the inferred structure and the reference 
solar model (model S) are quite small, generally less than 1%. The small 
differences provide a justification for the linearization procedure, based on 
the variational principle. This also means that the modern standard model 
of the Sun [101] provides an accurate description of the solar properties com- 
pared to the earlier solar model [106], used for the asymptotic inversions 
(Fig. 21). A significant improvement in the solar modeling was achieved by 
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Fig. 29. The relative differences between the Sun and the standard solar model 
[101]in the squared sound speed, c 2 , the adiabatic exponent, 7, the density, p, 
and the parameter of convective stability, A*, inferred from the solar frequencies 
determined from the 360-day series of SOHO MDI data. 



using more accurate radiative opacity data and by including the effects of 
gravitational settling of heavy elements and element diffusion. However, re- 
cent spectroscopic estimates of the heavy element abundance on the Sun, 
based on radiative hydrodynamics simulations of solar convection, indicated 
that the heavy element abundance on the Sun may be lower than the value 
used in the standard model [117]. The solar model with a low heavy element 
abundance do not agree with the helioseismology measurements (e.g. [118]). 
This problem in the solar modeling has not been resolved. Thus, the helio- 
seismic inferences of the solar structure lead to better understanding of the 
structure and evolution of the star, and have important applications in other 
fields of astrophysics. 

The prominent peak of the squared sound speed, Sc 2 /c 2 , at the base of 
the convection zone, r/R w 0.7, indicates on additional mixing which may 
be caused by rotational shear flows or by convective overshoot. The variation 
in the sound speed in the energy-generating core at r/R < 0.2 might be also 
caused by a partial mixing. 

The monotonic decrease of the adiabatic exponent, 7, in the core was 
recently explained by the relativistic corrections to the equation of state 
[119]. Near surface variations of 7, in the zones of ionization of helium and 
hydrogen, and below these zones, are most likely caused by deficiencies in the 
theoretical models of the weakly coupled plasma employed in the equation of 
state calculations [120]. 
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The monotonic decrease of the squared sound speed variation in the con- 
vection zone (r/R > 0.7) is partly due to an error in the solar seismic radius 
used to calibrate the standard model [107], and partly due to the inaccurate 
description of the subsurface layers by the standard solar model, based on a 
mixing-length convection theory. 



5.7 Regularized least-squares method 

The Regularized Least-Squares (RLS) method [116] is based on minimization 
of the quantity 



E 



i 



'n.l 



Uj(n,l) 



K (nd)8.f ( ni j) 5g 



dr, (125) 



in which the unknown structure correction functions, Sf / f and Sg/g, arc both 
represented by piece- wise linear functions or by cubic splines. The second 
integral specifies smoothness constraints for the unknown functions, in which 
Li and L 2 are linear differential operators, e.g. L 12 = d 2 /d 2 r; <Ji are error 
estimates of the relative frequency differences. 

In this inversion method, the estimates of the structure corrections are, 
once again, linear combinations of the frequency differences obtained from 
observations, and corresponding averaging kernels exist too. However, unlike 
the OLA kernels A(r$\ r), the RLS averaging kernels may have negative side- 
lobes and significant peaks near the surface, thus making interpretation of 
the inversion results to some extent ambiguous. Nevertheless, it works well 
in most cases, and may provide a higher resolution compared to the OLA 
method. 



5.8 Inversions for solar rotation 

The eigenfrequencies of a spherically-symmetrical static star are degenerate 
with respect to the azimuthal number m. Rotation breaks the symmetry 
and splits each mode of radial order, n, and angular degree, I, into (21 + 
1) components of m = —I,..., I (mode multiplets). The rotational frequency 
splitting can be computed using a more general variational principle derived 
by Lynden-Bcll and Ostriker [121]. From this variational principle, one can 
obtain mode frequencies uj n i m relative to the degenerate frequency to n i of the 
non-rotating star: 

Auj nlm = uj nlm -u n i = j- I [m£ ■ C + ie Q (i x £*)] QpdV, (126) 

+ nl JV 
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where ea is the unit vector defining the rotation axis, and Q = f2(r, 9) is the 
angular velocity which is a function of radius r and co-latitude 0, and I n i is 
the mode inertia. 



Equation (126) can be rewritten as a two-dimensional integral equation 
for Q(r,e): 



Here and r\ n l are the radial and horizontal components of cigenfunctions 
of the mean spherically symmetric structure of the Sun, P[ n (9) is an asso- 
ciated normalized Legendre function (./^(-P™) 2 sin OdO = 1). The kernels are 
symmetric relative to the equator, 9 — tt/2. Therefore, the frequency split- 
tings are sensitive only to the symmetric component of rotation in the first 
approximation. The non-symmetric component can, in principle, be deter- 
mined from the second-order correction to the frequency splitting, or from 
local helioseismic techniques, such as time-distance seismology. 

For a given set of observed frequency splitting, Auj n i m , eq. (126) con- 
stitutes a two-dimensional linear inverse problem for the angular velocity, 
f2(r, 9), which can be solved by the OLA or RLS techniques. 

5.9 Results for Solar Rotation 

As an example, we present the inversion results for solar rotation obtained 
from SOHO/MDI data. The frequency splitting data were obtained from the 
144-day MDI time series by J. Schou for j = 1, 36 and 1 < I < 250 [122]. 
The total number of measurements in this data set was M — 37366. 

Figure 30 shows results of inversion of the SOI-MDI data by the two 
methods. The results are generally in good agreement in most of the area 
where good averaging kernels were obtained. However, the results differ in 
the high-latitude region. In particular, a prominent feature of the RLS in- 
version at coordinates (0.2, 0.95) in Fig. 30a, which can be interpreted as a 
'polar jet', is barely visible in Fig. 30b, showing the OLA inversion of the 
same data. Therefore, obtaining reliable inversion results in this region and 
also in the shaded area is one of the main current goals of helioseismology. 
This can be achieved by obtaining more accurate measurements of rotational 
frequency splitting and improving inversion techniques. Of course, the rad- 
ical improvement can be made by observing the polar regions of the Sun. 




(127) 



where K^^r.9) represent the rotational splitting kernels: 





Fig. 30. Contour lines of the rotation rate (in nHz) inside the Sun obtained by 
inverting the rotational frequency splittings from a 144-day observing run from 
SOHO MDI by the RLS and SOLA methods. The shaded areas are the areas where 
the localized averaging kernels substantially deviate from the target positions. 



These measurements can be done by using spacecraft with an orbit highly 
inclined to the ecliptic plane, such as a proposed Solar Polar Imager (SPI) 
and POLARIS missions [123]. 

The most characteristic feature of solar rotation is the differential rota- 
tion of the convection zone, which occupies the our 30% of the solar radius. 
While the radiative core rotates almost uniformly, the equatorial regions of 
the convection zone rotate significantly faster than the polar regions. The 
main interest is in understanding the role of Sun's internal rotation in the 
dynamo process of generation of solar magnetic fields and the origin of the 11- 
year sunspot cycle. The results of these measurements (Fig. 31a) reveal two 
radial shear layers at the bottom of the convection zone (so-called tachocline) 
and in the upper convective boundary layer. A common assumption is that 
the solar dynamo operates in the tachocline area (interface dynamo) where 
it is easier to explain storage of magnetic flux than in the upper convec- 
tion zone because of the flux buoyancy. However, there are theoretical and 
observational difficulties with this concept. First, the magnetic field in the 
tachocline must be quite strong, ~ 60 — 160 kG, to sustain the action of 
the Coriolis force transporting the emerging flux tubes into high-latitude re- 
gions [124]. The magnetic energy of such field is above the equipartition level 
of the turbulent energy. Second, the back-reaction such strong field should 
suppress turbulent motions affecting the Reynolds stresses. Since these tur- 
bulent stresses support the differential rotation one should expect significant 
changes in the rotation rate in the tachocline. However, no significant varia- 
tions with the 11-year solar cycle are detected. Third, magnetic fields often 
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Fig. 31. The solar rotation rate as a function of radius at three latitudes. The 
horizontal lines indicate the rotation rate of the surface magnetic flux at the end of 
solar cycle 22 ("old magnetic flux") and at the beginning of cycle 23 ("new magnetic 
flux") [126] 



tend to emerge in compact regions on the solar surface during long periods 
lasting several solar rotations. This effect is known as "complexes of activity" 
or "active longitudes". However, the helioseismology observations show that 
the rotation rate of the solar tachocline is significantly lower than the surface 
rotation rate. Thus, magnetic flux emerging from the tachocline should be 
spread over longitudes (with new flux lagging the previously emerged flux) 
whether it remains connected to the dynamo region or disconnected. It is 
well-known that sunspots rotate faster than surrounding plasma. This means 
that the magnetic field of sunspots is anchored in subsurface layers. Obser- 
vations show that the rotation rate of magnetic flux matches the internal 
plasma rotation in the upper shear layer (Fig. 31b) indicating that this layer 
is playing an important role in the solar dynamo, and causing a shift in the 
dynamo paradigm [125]. 

Variations in solar rotation clearly related to the 11-year sunspot cycle are 
observed in the upper convection zone. These are so-called 'torsional oscilla- 
tions' which represent bands of slower and faster rotation, migrating towards 
the equator as the solar cycle progresses (Fig. 32). The torsional oscillations 
were first discovered on the Sun's surface [129], and then were found in the 
upper convection zone by helioseismology [130,131]. The depth of these evolv- 
ing zonal flows is not yet established. However, there are indications that they 
may be persistent through most of the convection zone, at least, at high lat- 
itudes [128]. The physical mechanism is not understood. Nevertheless, it is 
clear that these zonal flows are closely related to the internal dynamo mech- 
anism that produces toroidal magnetic field. On the solar surface, this field 
forms sunspots and active regions which tend to appear in the areas of shear 
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Fig. 32. a) Migration of the subsurface zonal flows with latitude during solar cycle 
23 from SOHO/MDI data [127]. Red shows zones of faster rotation, green and 
blue show slower rotation, b) Variations of the zonal flows with depth and latitude 
during the first 4 years after the solar minimum. [128] 



flows at the outer (relative to the equator) part of the faster bands. Thus, the 
torsional flows are an important key to understanding the solar dynamo, and 
one of the challenges is to establish their precise depth and detect correspond- 
ing variations in the thermodynamic structure of the convection zone. Recent 
modeling of the torsional oscillations by the Lorentz force feedback on differ- 
ential rotation showed that the poleward-propagating high-latitude branch 
of the torsional oscillations can be explained as a response of the coupled dif- 
ferential rotation/meridional flow system to periodic forcing in midlatitudes 
of either mechanical (Lorentz force) or thermal nature [132]. However, the 
main equatorward-propagating branches cannot be explained by the Lorenz 
force, but maybe driven by thermal perturbations caused by magnetic field 
[133]. It is intriguing that starting from 2002, during the solar maximum, the 
helioseismology observations show new branches of "torsional oscillations" 
migrating from about 45° latitude towards the equator (Fig. 32a). They in- 
dicate the start of the next solar cycle, number 24, in the interior, and are 
obviously related to magnetic processes inside the Sun. However, magnetic 
field of the new cycle appeared on the surface only in 2008. 
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6 Local-area helioseismology 

6.1 Basic principles 

In the previous sections we discussed methods of global helioseismology, which 
are based on inversions of accurately measured frequencies and frequency 
splitting of normal oscillation modes of the Sun. The frequencies are measured 
from long time series of observations of the Doppler velocity of the solar disk. 
These time series are much longer than the mode lifetimes, typically, two or 
three 36-day-long 'GONG months', that is 72 or 108 days. The long time 
series allow us to resolve individual mode peaks in the power spectrum, and 
accurately measure the frequencies and other parameters of these modes. 
However, because of the long integration times global helioseismology cannot 
capture the fast evolution of magnetic activity in subsurface layers of the 
Sun. Also, it provides only information about the axisymmetrical structure 
of the Sun and the differential rotation (zonal flows) . 

Local helioseismology attempts to determine the subsurface structure and 
dynamics of the Sun in local areas by analyzing local characteristics of solar 
oscillations, such as frequency and phase shifts and variations in wave travel 
times. This is a relatively new and rapidly growing field. It takes advantage 
of high-resolution observations of solar oscillations, currently available from 
the GONG+ helioseismology network and the space mission SOHO, and are 
anticipated from the SDO mission. 

6.2 Ring-diagram analysis 

Local helioseismology was pioneered by Douglas Gough and Juri Toomre [46] 
first proposed to measure oscillation frequencies of solar modes as a function 
of the wavevector, ui(k), (the dispersion relation) in local areas, and use these 
measurements for diagnostics of the local flows and thermodynamic proper- 
ties. They noticed that subsurface variations of temperature cause change in 
the frequencies, and that subsurface flows result in distortion of the dispersion 
relation because of the advection effect. 

This idea was implemented by Frank Hill [47] in the form of a ring- diagram 
analysis. The name of this technique comes from the ring appearance of the 
3D dispersion relation, uj — uj(k Xl ky), in the (k XJ ^y) pl&ns, where k x and hy 
are x- and y-components of the wave vector, k (Fig. 33). The ridges in the 
vertical cuts represent the same mode ridges as in Fig. 3, corresponding to 
the normal oscillation modes of different radial orders n. 

In the presence of a horizontal flow field, U = (U x ,U y ) the dispersion 
relation has the form: 

u = uj (k) + k-U = Lu a + (U x k x + Uyky), (129) 

where ooo(k) is the symmetrical part of the dispersion relation in the (k x ,k y )- 
plane. It depends only on the magnitude of the wave vector, fe. The power 
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Fig. 33. Three-dimensional power spectrum of solar oscillations, P(k x , k y ,uj). The 
vertical panels with blue background show the mode ridge structure similar to the 
global oscillation spectrum shown in Fig. 3. The horizontal cut with transparent 
background shows the ring structure of the power spectrum at a given frequency 
(courtesy of Amara Graps). 



spectrum, P(u>, k) for each k is fitted with a Lorentian profile [134]: 



P(u,k) 



A 



(uj-ujo + k x U x + k y Uy) 2 + T 2 fc 3 ' 



(130) 



where A,uo,r, and &o are respectively the amplitude, central frequency, line 
width and a background noise parameter. 

In some realizations, the fitting formula includes the line asymmetry 
(Sec. 3). Also, the central frequency can be fitted by assuming a power-law 
relation: ujq = ck p , where c and p are constants [135,47]. This relationship 
is valid for a polytropic adiabatic stratification, where p — 1/2 [46]. If the 
flow velocity changes with depth then the parameter, U, represent a velocity, 
averaged with the depth with a weighting factor proportional to the kinetic 
energy density of the waves, p£ ■ £ [136]: 



U 



(131) 



where £(z) = (£ r , £/j) is the wave amplitude, given by the mode displacement 
eigenfunctions (15. The integral is taken over the entire extent of the solar 
envelope. Equation (131) is solved by the RLS or OLA techniques (Sec. 5). 

The ring-diagram method has provided important results about the struc- 
ture and evolution of large-scale and meridional flows and dynamics of active 
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regions [137,134,138,127,139]. In particular, large-scale patterns of subsurface 
flows converging around magnetic active regions were discovered [138]. These 
flows cause variations of the mean meridional circulation with the solar cycle 
[134], which may affect transport of magnetic flux of decaying active regions 
from low latitudes to the polar regions, and thus change the duration and 
magnitude of the solar cycles [140]. 

However, the ring-diagram technique in the present formulation has lim- 
itations in terms of the spatial and temporal resolution and the depth cov- 
erage. The local oscillation power spectra are typically calculated for regions 
with the horizontal size covering 15 heliographic degrees (~ 180 Mm). This 
is significantly larger than the typical size of supergranulation and active 
regions (~ 30 Mm). There have been attempts to increase the resolution 
by doing the measurements in overlapping regions (so-called "dense-packed 
diagrams"). However, since such measurements are not independent, their 
resolution is unclear. The measurements of the power spectra calculated for 
smaller regions (2-4 degrees in size) increase the spatial resolution but de- 
crease the depth coverage [141]. 

6.3 Time-distance helioseismology (Solar tomography) 

Further developments of local seismology led to the idea to perform mea- 
surements of local wave distortions in the time-distance space instead of the 
traditional frequency- wavenumber Fourier space [48]. In this case, the wave 
distortions can be measured as perturbations of wave travel times. However, 
because of the stochastic nature of solar waves it is impossible to track indi- 
vidual wave fronts. Instead, it was suggested to use a cross-covariance (time- 
distance) function that provides a statistical measure of the wave distortion. 
Indeed, by cross-correlating solar oscillation signals at two points one may 
expect that the main contribution to this cross-correlation will be from the 
waves traveling between these points along the acoustic ray paths [142,143]. 
Thus, the cross-ccovariance function calculated for oscillation signals mea- 
sured at two points separated by a distance, A, for various time lags, r, has 
a peak when the time lag is equal to the travel time of acoustic waves be- 
tween these points. Physically, the cross-covariance function corresponds to 
the Green's function of the wave equation, representing the wave signal from 
a point source. Of course, in reality, because of the finite wavelength effects, 
non-uniform distribution of acoustic sources, and complicated wave interac- 
tion with turbulence and magnetic fields the interpretation of the travel-time 
measurements is extremely challenging. Various approximations are used to 
relate the observed perturbations of the travel times to the internal proper- 
ties such as sound-speed perturbations and flow velocities. We discuss the 
basic principles and the current status of the time-distance helioseismology 
method in Sec. 7. 



60 Alexander G. Kosovichev 



6.4 Acoustic holography and imaging 

The acoustic holography [144] and acoustic imaging [51] techniques are de- 
veloped on the principles of day-light imaging by collecting over large areas 
on the solar surface coherent acoustic signals emitted from selected target 
points of the interior. The idea is that the constructed this way signals con- 
tain information about objects located below the surface because of wave 
absorption or scattering at the target points. The phases of individual sig- 
nals are calculated by using the time-distance relation, t(A), f or acoustic 
waves traveling along the ray paths. The constructed signals, VwunW, are 
calculated using the following relation [145]: 

T 2 

where tp(A,t + r) is the azimuthal-avcraged signal at a distance A from a 
target point at time t±r(A). The summation variable r is equally spaced in 
the interval (ri,T2); and the weighting factor, W oc (sin A/t 2 ) 1 / 2 , describes 
the geometrical spreading of acoustic waves with distance. The positive sign in 
equation (132) corresponds to Vw constructed with waves traveling outward 
from a target point ("egression signal" [144]), while the negative sign provides 
ipi n constructed with the incoming waves ("ingression signal"). 

The amplitude and phase of the constructed signals contain information 
about subsurface perturbation. A practical approach to extract this is to 
cross-correlate the outgoing and incoming signals [146,147]: 

C(t) = J Vi„(i')V'out(i' + t)dt', (133) 

and then to measure time shifts of this function for various target positions 
relative to the corresponding quiet Sun values. These measurements corre- 
spond to the travel-time variations obtained by time-distance helioseismol- 
ogy [148,149]. Further analysis of the travel-time variations is similar to the 
time-distance helioseismology method [50]. The advantages and disadvan- 
tages of the time-distance helioseismology and acoustic holography/imaging 
are not clear. Both, approaches are being tested using various types of artifi- 
cial data and applied for measuring subsurface structures and flows. Most of 
the current inferences of subsurface structures and flows have been obtained 
using the time-distance approach [48,50]. The time-distance helioseismology 
method, also called solar tomography is described in more detail in the fol- 
lowing section. 

7 Solar tomography 
7.1 Time-distance diagram 

Solar acoustic waves (p-modes) are excited by turbulent convection near the 
solar surface and travel through the interior with the speed of sound. Because 
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the sound speed increases with depth the waves are refracted and reappear 
on the surface at some distance from the source. The wave propagation is 
illustrated in Figure 34. Waves excited at point A will reappear at the sur- 
face points B, C, D, E, F, and others after propagating along the ray paths 
indicated by the curves connecting these points. 



E I 

Fig. 34. A cross-section diagram through the solar interior showing a sample of 
wave paths inside the Sun. 

The basic idea of time-distance helioseismology, or helioseismic tomogra- 
phy, is to measure the acoustic travel time between different points on the 
solar surface, and then to use these measurements for inferring variations of 
wave-speed perturbations and flow velocities in the interior by inversion [48] . 
This idea is similar to seismology of Earth. However, unlike in Earth, the 
solar waves are generated stochastically by numerous acoustic sources in a 
subsurface layer of turbulent convection. 

Therefore, the wave travel time is determined from the cross-covariance 
function, \P(t, A), of the oscillation signal, f{t,r): 



where A is the horizontal distance between two points with coordinates r\ 
and r-i, t is the lag time, and T is the total time of the observations. The 
normalized cross-covariance function is called cross-correlation. The time- 
distance analysis is based on non-normalized cross-covariance. Because of 
the stochastic nature of solar oscillations, function W must be averaged over 
some areas to achieve a good signal-to-noise ratio sufficient for measuring 
the travel times. The oscillation signal, f(t, r), is measured from the Doppler 
shift or intensity of a spectral line. A typical cross-covariance function ob- 
tained from full-disk solar observations of the Doppler shift shown in Fig. 35a 
displays a set of ridges. The ridges correspond to acoustic wave packets trav- 



A 





(134) 
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cling between two points on the surface directly through the interior or with 
intermediate reflections (bounces) from the surface as illustrated in Figure 
34 



Observational Time-Distance Diagram Theoretical Time-Distance Diagram 




50 100 150 50 100 150 

distance, deg distance, deg 

a) b) 



Fig. 35. The observational (a) and theoretical (b) cross-covariance functions (time- 
distance diagrams) as a function of distance on the solar surface, A, and the delay 
time, r. The lowest set of ridges ('first bounce') corresponds to waves propagated to 
the distance, A, without additional reflections from the solar surface. The second 
from the bottom ridge ('second bounce') is produced by the waves arriving to the 
same distance after one reflection from the surface, and the third ridge ('third 
bounce') results from the waves arriving after two bounces from the surface. The 
backward ridge at r « 250 min is a continuation of the second-bounce ridge due 
to the choice of the angular distance range from to f80 degrees (that is, the 
counterclockwise distance ADF in Fig. 34 is substituted with the clockwise distance 
AF). Because of foreshortening close to the solar limb the observational cross- 
covariance function covers only ~ f fO degrees of distance. 



The waves originated at point A may reach point B directly (solid curve) 
forming the first-bounce ridge, or after one bounce at point C (dashed curve) 
forming the second-bounce ridge, or after two bounces (dotted curve) - the 
third-bounce ridge and so on. Because the sound speed is higher in the deeper 
layers the direct waves arrive first, followed by the second-bounce and higher- 
bounce waves. 

The cross-covariance function represents a time- distance diagram, or a 
solar 'seismogram'. Figure 36 shows the cross-covariance signal as a function 
of time for the travel distance, A, of 30 degrees. It consists of three wave 
packets corresponding to the first, second and third bounces. Ideally like in 
Earth seismology, the seismogram can be inverted to infer the structure and 
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flows using a wave theory. However, in practice, modeling the wave fronts 
is a computationally intensive task. Therefore, the analysis is performed by 
measuring and inverting the phase and group travel times of the wave packets 
employing various approximations, the most simple and powerful of which is 
the ray-path approximation. 



1.0 




50 100 150 200 

time, min 



Fig. 36. The observed cross-covariance signal as a function of time at the distance 
of 30 degrees. 



Generally, the observed solar oscillation signal corresponds to displace- 
ment or pressure perturbation, and can be represented in terms of the nor- 
mal modes eigenfunctions. Therefore, the cross-covariance function also can 
be expressed in terms of the normal modes. In addition, it can be represented 
as a superposition of traveling wave packets, as we show in the next subsec- 
tion [50] . An example of the theoretical cross-covariance function calculated 
using normal p-modes of the standard solar model is shown in Fig. 35b. This 
model reproduces the observational cross-covariance function very well in the 
observed range of distances, from to 90 degrees. The theoretical model was 
calculated for larger distances than the corresponding observational diagram 
in Fig. 35a, including points on the far side of the Sun, which is not acces- 
sible for measurements. A backward propagating ridge originating from the 
second-bounce ridge at 180 degrees is a geometrical effect due to the choice 
of the range of the angular distance from to 180 degrees. In the theoreti- 
cal diagram (Fig. 35b) one can notice a very weak backward ridge between 
30 and 70 degrees and at 120 min. This ridge is due to reflection from the 
boundary between the convection and radiative zones. However, this signal 
has not been detected in observations. 
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7.2 Wave travel times 

For simplicity we consider solar oscillation signals observed not far from the 
disk center and describe these in terms of the radial displacement neglecting 
the horizontal displacement. The general theory was developed by Nigam 
and Kosovichev [150]. In the simple case, the solar oscillation signal can be 
represented in terms of the radial eigenfunctions (42): 

f(t, r, 6,<f>) = J2 am m & n ' l ' m) (r, 9, <j>) cxp(iuj nlm t + t<t> nlm ), (135) 

nlm 

where n, I and m are the radial order, angular degree and angular order of a 
normal mode respectively, ^„; m (r, 9, <p) is a mode eigenfunction in the spher- 
ical coordinates, r, 9 and (j>, uj n im is the eigenfrequency, and (\> n im is an initial 
phase of the mode. Using equation (135), we calculate the cross-covariancc 
function, and express it as a superposition of traveling wave packets. Such 
a representation is important for interpretation of the time-distance data. A 
similar correspondence between the normal modes and the wave packets has 
been discussed for surface oscillations in Earth's seismology [151] and also 
for ocean waves [152]. 

To simplify the analysis, we consider the spherically symmetrical case. In 
this case, the mode eigenfrequencies do not depend on the azimuthal order 
m. For a radially stratified sphere, the eigenfunctions can be represented in 
terms of spherical harmonics Yi m (9, cf)) (42): 

£< n ''' m) (r,M) =^ n ' l \r)Y lm {e,<j>), (136) 

where £r"'^( r ) i s th e radial eigenfunction [153]. 

Using, the convolution theorem [154] we express the cross-covariance func- 
tion in terms of a Fourier intergral: 

/oo 
F(oj,r 1 )F*(w,r 2 )exp(iojT)duj, (137) 

where F(ui,r) is Fourier transform of the oscillation signal f(t,r). 

The oscillation signal is considered as band-limited and filtered to select 
a p-mode frequency range using a Gaussian transfer function: 



G(uj) = exp 



1 / uj — UJq 



2 V Suj 



(138) 



where u> is the cyclic frequency, tu is the central frequency and Suj is the 
characteristic bandwidth of the filter. The cross-covariance function in Fig. 1 
displays three sets of ridges which correspond to the first, second and third 
bounces of packets of acoustic wave packets from the surface. 
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The time series used in our analysis are considerably longer than the 
travel time r, therefore, we can neglect the effect of the window function, 
and represent F(ui,r) in the form 



%r,M)»^ [r)Y lm (6, cj>)5{u - u nl ) exp 



nlm 



1 / Lt) — OJq 

2 { 5uj 



(139) 

where 5{x) is the delta- function, u n i are frequencies of the normal modes, and 
A is the amplitude of the Gaussian envelope of the amplitude spectrum at 
uj = luq- In addition, we assume the normalization conditions: £,i n ' l \n) = 1, 
a n i = AG{uj). Then, the cross-covariance function is 



<P(r, A) = A 2 ^2 ex P 



nl 



Slu 



lUJnlT 



m— — l 



(140) 

where 6\ , <j>\ and &i , <p2 are the spherical heliographic coordinates of the two 
observational points. The sum of the spherical function products 



Yim(0i,<tn)Y£ n (e2,(h)=aiPi{co8A), 



(141) 



m— — I 



where P;(cosZ\) is the Legendre polynomial, A is the angular distance be- 
tween points 1 and 2 along the great circle on the sphere, cos A = cos 6\ cos #2+ 
sin 0i sin 2 cos(</>2 — </>i), and on = ^ Air/ (21 + 1). Then, the cross-covariance 
function is: 



&(t, A)^ A 2 Y^ ottPi (cos A) exp 

nl 

For large values of I A, but when A is small 
Pi (cos A) ~ -\ / ^ ; cos 

Thus, 



(5w 



(142) 



(143) 



!^(r,Zi)=A 2 ^-^=exp 



(w n ; - w ) 



2 1 



<Su; 2 



cos(w n iT) cos(LZ\). (144) 



Now the double sum can be reduced to a convenient sum of integrals if we 
regroup the modes so that the outer sum is over the ratio v = uj n i / L and the 
inner sum is over w„; . 

According to the ray-path theory, the travel distance A of an acoustic 
wave is determined by the ratio v, which represent the horizontal angular 
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phase velocity (v = cj n i/L = (co n i/kh)/r). Because of the band-limited nature 
of the function G, only values of L which are close to Lq = luq/v contribute 
to the sum. We consider the relation L vs u n i as a continuous function along 
the mode ridges (Fig. 3), and expand L near the central frequency ojq: 



L - L + a K^nl - W ) = 1 

OUJnl V U 



W 



where u = duj n i/dL. Furthermore, 
cos(w„;)r) cos(LZi) = COS 



A 



T — 



Aojq 



(145) 



(146) 



and the other term is identical except that t has been replaced with — r 
(negative time lag). The result is that the double sum in equation (146) 
becomes 



*(t,A)^A 2 J2 



l Va 



J^exp 



cos 



±T - 



A 



+ 



- - - I Awo 



(147) 

The inner sum can be approximated by an integral, considering ui n i as a 
continuous variable along the mode ridges: 



/oo 
dui exp 
-oo 

V 7r Sou 2 exp 



(u - UJq) 2 
5^ 



cos 
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u 



5uj 2 














cos 


UJQ 


H)] 


4 









Auj a = 
(148) 



The integration limits reflect the fact that the amplitude function G{uS) is 
essentially zero for very large and very small frequencies. Finally, the cross- 
covariance is expressed in the following form [50] : 



\P(t, A) = B ^2 cos I w o (t - r p h)] exp 



Suj 2 



'grf 



(149) 



where B is constant, r p h — A/v and r gr = A/u are the phase and group 
travel times. Equation (149) has the form of a Gabor wavelet. The phase and 
group travel times are measured by fitting individual terms of equation (149) 
to the observed cross-covariance function using a least-squares technique. 



7.3 Deep- and surface-focus measurement schemes 

As we have pointed out the travel-time measurements require averaging of 
the cross-covariance function in order to obtain a good signal-to-noise ratio. 
Two typical schemes of the spatial averaging suggested by Duvall [155] arc 
shown in Fig. 37. 
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Surface focusing 




Fig. 37. The regions of ray propagation (shaded areas) as a function of depth, z, 
and the radial distance, A, from a point on the surface for two observing schemes: 
'surface focusing' (a) and 'deep focusing' (b). The rays are also averaged over a 
circular regions on the surface, forming three-dimensional figures of revolution. 



For the so-called 'surface-focusing' scheme (Fig. 37a) the measured travel 
times are mostly sensitive to the near surface condition at the central point 
where the ray paths are focused. However, by measuring the travel times for 
several distances and applying an inversion procedure it is possible to infer 
the distribution of the variations of the wave speed and flow velocities with 
depth. The averaging also can be done in such a way that the 'focus' point 
is located beneath the surface. An example of the 'deep-focusing' scheme is 
shown in Fig. 37b. In this case the travel times are more sensitive to deep 
structures but still inversions are required for correct interpretation. 

7.4 Sensitivity kernels: Ray-path approximation 

The travel-time inversion procedures are based on theoretical relations be- 
tween the travel-time variations and interior properties constituting the for- 
ward problem of local helioseismology. Similarly to global helioseismology, 
these relations are expressed in the form linear integral equations with sen- 
sitivity kernels. Two basic types of the sensitivity kernels have been used: 
ray-path kernels [50] and Born-approximation kernels [156-158]. The ray- 
path kernels are based on a simple and generally robust theoretical ray ap- 
proximation, but they do not take into account finite wavelength effects and 
thus are not sufficiently accurate for diagnostics of small-scale structures. For 
reliable inferences it is important to use both these kernels. 

In the ray approximation, the travel times are sensitive only to the per- 
turbations along the ray paths given by Hamilton's equations (72). The vari- 
ations of the phase travel time obey the Fermat's Principle: 




(150) 
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where 5k is the perturbation of the wave vector, fc, due to the structural 
inhomogeneities and flows along the unperturbed ray path, r. Using the dis- 
persion relation for acoustic waves in the convection zone the travel-time vari- 
ations can be expressed in terms of the sound-speed, magnetic field strength 
and flow velocity. 

The dispersion relation for magnetoacoustic waves in the convection zone 

is 

(w - k ■ Uf = to 2 c + k 2 c 2 fl (151) 
where U is the flow velocity, uj c is the acoustic cut-off frequency, c 2 — 



\[c 2 + c\ + 



\J (c 2 + c 2 A ) 2 — 4c 2 (fc • ca) 2 /k 2 ^j is the fast magnetoacoustic speed, 



ca = B/^/AiFp is the vector Alfven velocity, B is the magnetic field strength, 
c is the adiabatic sound speed, and p is the plasma density. If we assume that, 
in the unperturbed state U = B = 0, then, to the first-order approximation 



St 



f Un-U) | 5c s | / Suj c \ ul | 1 (c\ (k-c A ) 
Jr [ c 2 c \u c ) lj 2 c 2 S 2 V 



ds, 



k 2 c 2 

" (152) 

where n is a unit vector tangent to the ray, S = k/u; is the phase slowness. 

Then, we separate the effects of flows and structural perturbations by 
measuring the travel times of acoustic waves traveling in opposite directions 
along the same ray path, and calculating the difference, Tdiff and the mean, 
T m ean, of these reciprocal travel times: 



Sms = -2 



/ { ^ds; (153) 
Jr c 



St„ 



5c s ( SujA uj 2 + 1{ca_ (fc ■ ca, | g 



c \ oj r I lo 2 c 2 S 2 V c 2 k 2 c 2 



ds. (154) 



Anisotropy of the last term of equation (154) allows us to separate, at least 
partly, the magnetic effects from the variations of the sound speed and the 
acoustic cut-off frequency. The acoustic cut-off frequency, uj c may be per- 
turbed by surface magnetic fields and by temperature and density inhomo- 
geneities. The effect of the cut-off frequency variation depends strongly on 
the wave frequency, and, therefore, it results in a frequency dependence in 

7"mcan • 

In practice, the travel times are measured for from the cross-covariance 
functions between selected central points on the solar surface and surround- 
ing quadrants symmetrical relative to the North, South, East and West di- 
rections. In each quadrant, the travel times are averaged over narrow ranges 
of the travel distance, A The travel times of the northward-directed waves 
are subtracted from the times of the south-directed waves to yield the time, 
r diff' which is predominantly sensitive to subsurface north-south flows. Sim- 
ilarly, the time differences, t™ , between westward- and eastward directed 
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waves yields a measure of the east- ward flows. The time, t^[ s , between the 
outward- and inward-directed waves, averaged over the full annuli, is mainly 
sensitive to vertical flows and divergence of the horizontal flows. This repre- 
sents a cross-talk effect between the vertical flows and horizontal flows, which 
is difficult to resolve when the vertical flows are weak [159]. 

Thus, the effects of flows and structural perturbations are separated from 
each other by taking the difference and the mean of the reciprocal travel 
times: 

Sr dis « -2 / { -^fids; (155) 
Jr c 



[ —Sds, (156) 
Jr c 



where c is the adiabatic sound speed, n is a unit vector tangent to the ray, 
S — k/uj is the phase slowness, 5w is the local wave speed perturbation: 

Sw _5c , (c 2 A (kc A ) 2 \ 

c ~ c + Hc 2 J' (i57j 

Magnetic field causes anisotropy of the mean travel times, which allows us to 
separate, in principle, the magnetic effects from the variations of the sound 
speed (or temperature). So far, only a combined effect of the magnetic fields 
and temperature variations has been measured reliably. 



7.5 Born approximation 

The development of a more accurate theory for the travel times, based on 
the Born approximation is currently under way [156,160,161,157,158]. 

One unexpected feature of the single-source travel-time kernels calculated 
in the Born approximation is that these kernels have zero value along the ray 
path (called 'banana-doughnut kernels'). Examples of the Born kernels for 
the first and the second bounces are shown in Fig. 38. The kernels are mostly 
sensitive to perturbations within the first Fresnel zone. 

Figure 39 shows the test results for both the ray and Born approximations 
for a simple model of a smooth sphere in an uniform medium by comparing 
with precise numerical results [160]. These results show that for typical per- 
turbations in the solar interior the Born approximation is sufficiently accu- 
rate, while the ray approximation significantly overestimates the travel times 
for perturbations smaller than the size of the first Fresnel zone. That means 
that the inversion results based on the ray theory may underestimate the 
strength of the small-scale perturbations. The comparison of the inversion 
results for sub-surface sound-speed structures beneath sunspots have showed 
a very good agreement between the ray-paths and Born theories [158]. 
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A = 60 deg, second bounce 




Fig. 38. Travel-time sensitivity kernels in the first Born approximation for sound- 
speed variations as a function of the horizontal, x, and vertical, y, coordinates for: 
a) the first-bounce signal for distance A = 6 degrees, b) the second-bounce signal 
for A — 60 degrees. The solid curves show the corresponding ray paths at frequency 
v = 3 mHz [162]. 



8 Inversion results of solar acoustic tomography 

The results of test inversions (e.g. [50,159,55]) demonstrate an accurate re- 
construction of sound-speed variations and the horizontal components of sub- 
surface flows. However, vertical flows in deep layers are not resolved because 
of the predominantly horizontal propagation of the rays in these layers. The 
vertical velocities are also systematically underestimated in the upper layers. 
When the vertical flow is weak, e.g. such as in supergranulation, the vertical 
velocity is not estimated correctly, because the trave-time signal is dominated 
by the horizontal flow divergence. In such situation, it is difficult to deter- 
mine even the direction of the vertical flow [55]. Similarly, the sound-speed 
variations are underestimated in the deep layers and close to the surface. 
These limitations of the solar tomography should be taken into account in 
interpretation of the inversion results. 

Here, I briefly present some examples of the local helioseismology infer- 
ences obtained by inversion of acoustic travel times. 

8.1 Diagnostics of supergranulation. 

The data used were for 8.5 hours on 27 January, 1996 from the high resolution 
mode of the MDI instrument. The results of inversion of these data are shown 
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Fig. 39. Tests of the ray and Born approximations: travel times for smooth spheres 
as functions of sphere radius at half maximum. The solid lines are the numeri- 
cal results. The dashed curves are the Born approximation travel times and the 
dotted lines are the first order ray approximation. The left panel shows the two 
perturbations of the relative amplitude, A = ±0.05. The right panel is for the cases 
A = ±0.1. [160] 



in Figure 40 [50]. It has been found that, in the upper layers, 2-3 Mm deep, 
the horizontal flow is organized in supergranular cells, with outflows from 
the center of the supergranules. The characteristic size of the cells is 20- 
30 Mm. Comparing with MDI magnetograms, it was found that the cell 
boundaries coincide with the areas of enhanced magnetic field. These results 
are consistent with the observations of supcrgranulation on the solar surface. 
However, in the layers deeper than ~ 5 Mm, the supergranulation pattern 
disappears. The inversions show an evidence of reverse converging flows at 
the depth of <~ 10 Mm [159]. This means that supergranulation is a relatively 
shallow phenomenon. 

8.2 Structure and dynamics of sunspot 

The high-resolution data from the SOHO and Hinode space missions have al- 
lowed us to investigate the structure and dynamics beneath sunspots. Figure 
41 shows an example of the internal structure of a large sunspot observed on 
June 17, 1998 [163]. An image of the spot taken in the continuum is shown at 
the top. The wave-speed perturbations under the sunspot are much stronger 
than these of the emerging flux, and can reach <~ 3 km/s. It is interesting that 
beneath the spot the perturbation is negative in the subsurface layers and 
becomes positive in the deeper interior. One can suggest that the negative 
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Fig. 40. The supergranulation horizontal flow velocity field (arrows) and the sound- 
speed perturbation (color background) at the depths of 1.4 Mm (a) and 5.0 Mm 
(b), as inferred from the SOHO/MDI high-resolution data of 27 January 1996. [50] 



perturbations beneath the spot are, probably, due to the lower temperature. 
It follows that magnetic inhibition of convection that makes sunspots cooler 
is most effective within the top 2-3 Mm of the convection zone. The strong 
positive perturbation below suggests that the deep sunspot structure is hot- 
ter than the surrounding plasma. However, the effects of temperature and 
magnetic field have not been separated in these inversions. Separating these 
effects is an important problem of solar tomography. These data also show 
at a depth of ~ 4 Mm connections to the spot of small pores, which have the 
same magnetic polarity as the main spot. The pores of the opposite polarity 
are not connected to the main sunspot. This suggests that sunspots represent 
a tree-like structure in the upper convection zone. 

Figure 42 shows the subsurface structures and flows beneath a sunspot 
obtained from Hinode [164]. A vertical cut along the East- West direction 
approximately in the middle of a large sunspot observed in AR 10953, May 
2, 2007, (Fig. 42a), shows that the wave speed anomalies extend about half 
of the sunspot size beyond the sunspot penumbra into the plage area. In the 
vertical direction, the negative wave speed perturbation extends to a depth 
of 3-4 Mm. The positive perturbation is about 9 Mm deep, but it is not clear 
whether it extends further, because our inversion cannot reach deeper layers 
because of the small field of view. Similar two-layer sunspot structures were 
observed before from SOHO/MDI [163] (Fig. 41). But, it is striking that the 
new images strongly indicate on the cluster structure of the sunspot [165]. 
This was not previously seen in the tomographic images of sunspots obtained 
with lower resolution. 

The high-resolution flow field below the sunspot is also significantly more 
complicated than the previously inferred from SOHO/MDI [166], but reveals 
the same general converging downdraft pattern. A vertical view of an aver- 
aged flow field (Fig. 42b) shows nicely the flow structure beneath the active 
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Fig. 41. The sound-speed perturbation in a large sunspot observed on June 20, 
1998, are shown as vertical and horizontal cuts. The horizontal size of the box is 13 
degrees (158 Mm), the depth is 24 Mm. The positive variations of the sound speed 
are shown in red, and the negative variations (just beneath the sunspot)are in blue. 
The upper semitransparent panel is the surface intensity image (dark color shows 
umbra, and light color shows penumbra). In panel b) the horizontal sound-speed 
plane is located at the depth of 4 Mm, and shows long narrow structures ('fingers') 
connecting the main sunspot structure with surrounding pores the same magnetic 
polarity as the spot [163]. 



region. Strong downdrafts are seen immediately below the sunspot's surface, 
and extends up to 6 Mm in depth. A little beyond the sunspot's boundary, 
one can find both upward and inward flows. Clearly, large-scale mass circu- 
lations form outside the sunspot, bringing plasma down along the sunspot's 
boundary, and back to the photosphere within about twice of the sunspot's 
radius. It is remarkable that such an apparent mass circulation is obtained 
directly from the helioseismic inversions without using any additional con- 
straints, such as forced mass conservation. Previously, the circulation pattern 
was not that clear. 

8.3 Large-scale and meridional flows 

Time-distance helioseismology [167] and also local measurements of the p- 
mode frequency shifts by the 'ring-diagram' analysis [134,137,138], have pro- 
vided synoptic maps of subsurface flows over the whole surface of the Sun. 
Figure 43 shows a portion of a high-resolution synoptic flow map at the 
depth of 2 Mm below the surface. In addition, to the supergranulation pat- 
tern these maps reveal large-scale converging plasma flow around the active 
regions where magnetic field is concentrated. These flows are particularly 
well visible in low- resolution synoptic flow maps (Fig. 44). The characteristic 
speed of these flows is about 50 m/s. 

These stable long-living flow patterns affect the global circulation in the 
Sun. It is particularly important that these flows change the mean meridional 
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Fig. 42. Wave speed perturbation and flow velocities beneath sunspots from Hinode 
data[164] 

flow from the equator to the poles, slowing it down during the solar max- 
imum years (Fig. 45). This may have important consequences for the solar 
dynamo theories which invoke the meridional flow to explain the magnetic 
flux transport into the polar regions and the polar magnetic field polarity 
reversals usually happening during the period of maximum of solar activity. 

9 Conclusion and outlook 

During the past decade thanks to the long-term continuous observations from 
the ground and space the physics of solar oscillations made a tremendous 
progress in understanding the mechanism of solar oscillations, and in devel- 
oping new techniques for helioseismic diagnostics of the solar structure and 
dynamics. However, many problems are still unresolved. Most of them are 
related to phenomena in strong magnetic field regions and in the deep in- 
terior. The prime helioseismology tasks are to detect processes of magnetic 
field generation and transport in the solar interior, and formation of active 
regions and sunspots. This will be help to understand the physics of the solar 
dynamo and the cyclic behavior of solar activity. 
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Fig. 43. A portion of a synoptic subsurface flow map at depth of 2 Mm. The color 
background shows the distribution of magnetic field on the surface [167]. 
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Fig. 44. Subsurface synoptic flow maps at three depths. The color background 
shows the distribution of magnetic field on the surface [167]. 



For solving these tasks it is very important to continue developing realistic 
MHD simulations of solar convection and oscillations and to obtain contin- 
uous high-resolution helioseismology data for the whole Sun. The recent ob- 
servations from Hinode have convincingly demonstrated advantages of high- 
resolution helioseismology, but unfortunately such data are available only for 
small regions and for short periods of time. A new substantial progress in 
observations of solar oscillations is expected from the Solar Dynamics Obser- 
vatory (SDO) space mission launched in February 2010. 
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Fig. 45. Evolution of subsurface meridional flow during 1996-2002 for various Car- 
rington rotations [167]. 



The Helioseismic and Magnetic Imager (HMI) instrument on SDO will 
provide uninterrupted Doppler shift measurements over the whole visible disk 
of the Sun with a spatial resolution of 0.5 arcsec per pixel (4096 x 4096 images) 
and 40-50 sec time cadence. The total amount of data from this instrument 
will reach 2 Tb per day. This tremendous amount of data will be processed 
through a specially developed data analysis pipeline and will provide high- 
resolution maps of subsurface flows and sound-speed structures [53] . These 
data will enable investigations of the multi-scale dynamics and magnetism of 
the Sun and also contribute to our understanding of the Sun as a star. 

The tools that will be used in the HMI program include: helioseismology 
to map and probe the solar convection zone where a magnetic dynamo likely 
generates this diverse range of activity; measurements of the photosphcric 
magnetic field which results from the internal processes and drives the pro- 
cesses in the atmosphere; and brightness measurements which can reveal the 
relationship between magnetic and convective processes and solar irradiance 
variability. 
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Hclioseismology, which uses solar oscillations to probe flows and struc- 
tures in the solar interior, is providing remarkable new perspectives about 
the complex interactions between highly turbulent convection, rotation and 
magnetism. It has revealed a region of intense rotational shear at the base of 
the convection zone, called the tachocline, which is the likely seat of the global 
dynamo. Convective flows also have a crucial role in advecting and shearing 
the magnetic fields, twisting the emerging flux tubes and displacing the pho- 
tospheric footpoints of magnetic structures present in the corona. Flows of 
all spatial scales influence the evolution of the magnetic fields, including how 
the fields generated near the base of the convection zone rise and emerge at 
the solar surface, and how the magnetic fields already present at the surface 
are advected and redistributed. Both of these mechanisms contribute to the 
establishment of magnetic field configurations that may become unstable and 
lead to eruptions that affect the near-Earth environment. 

New methods of local-area helioseismology have begun to reveal the great 
complexity of rapidly evolving 3-D magnetic structures and flows in the sub- 
surface shear layer in which the sunspots and active regions are embedded. 
Most of these new techniques were developed during analysis of MDI observa- 
tions. As useful as they are, the limitations of MDI telemetry availability and 
the limited field of view at high resolution has prevented the full exploitation 
of the methods to answer the important questions about the origins of solar 
variability. By using these techniques on continuous, full-disk, high-resolution 
observations, HMI will enable detailed probing of dynamics and magnetism 
within the near-surface shear layer, and provide sensitive measures of varia- 
tions in the tachocline. 

The scientific operation modes and data products can be divided into four 
main areas: global helioseismology, local-area helioseismology, line-of-sight 
and vector magnetography and continuum intensity studies. The principal 
data flows and products are summarized in Figure 46. 

Global Helioseismology:Diagnostics of global changes inside the Sun. The 
traditional normal-mode method described in Sec. 4-5, will provide large-scale 
axisymmetrical distributions of sound speed, density, adiabatic exponent and 
flow velocities through the whole solar interior from the energy-generating 
core to the near-surface convective boundary layer. These diagnostics will 
be based on frequencies and frequency splitting of modes of angular degree 
up to 1000, obtained for several day intervals each month and up to Z=300 
for each 2-month interval. These will be used to produce a regular sequence 
of internal rotation and sound-speed inversions to allow observation of the 
tachocline and average near surface shear. 

Local- Area Helioseismology: 3D imaging of the solar interior. The new meth- 
ods of local-area helioseismology (Sec. 6-reftomography) , time-distance tech- 
nique, ring-diagram analysis and acoustic holography represent powerful tools 
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Fig. 46. A schematic illustration of the Solar Dynamics Observatory HMI data 
analysis pipeline and data products. The dark shaded area indicates Level- 1 data 
products. The boxes to the right of this area represent intermediate and final Level- 
2 data products. The data products are described in detail in the HMI Science Plan 
[53]. 



for investigating physical processes inside the Sun. These methods on mea- 
suring local properties of acoustic and surface gravity waves, such as travel 
times, frequency and phase shifts. They will provide images of internal struc- 
tures and flows on various spatial and temporal scales and depth resolution. 
The targeted high-level regular data products include: 

• Full-disk velocity and sound-speed maps of the upper convection zone 
(covering the top 30 Mm) obtained every 8 hours with the time-distance 
methods on a Carrington grid; 

• Synoptic maps of mass flows and sound-speed perturbations in the upper 
convection zone for each Carrington rotation with a 2-degree resolution, 
from averages of full disk time-distance maps; 

• Synoptic maps of horizontal flows in upper convection zone for each Car- 
rington rotation with a 5 degree resolution from ring-diagram analyses. 

• Higher-resolution maps zoomed on particular active regions, sunspots 
and other targets, obtained with 4-8-hour resolution for up to 9 days 
continuously, from the time-distance method; 
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• Deep-focus maps covering the whole convection zone depth, 0-200 Mm, 
with 10-15 degree resolution; 

• Far-side images of the sound-speed perturbations associated with large 
active regions every 24 hours. 

The HMI science investigation addresses the fundamental problems of 
solar variability with studies in all interlinked time and space domains, in- 
cluding global scale, active regions, small scale, and coronal connections. One 
of the prime objectives of the Living With a Star program is to understand 
how well predictions of evolving space weather variability can be made. The 
HMI investigation will examine these questions in parallel with the funda- 
mental science questions of how the Sun varies and how that variability drives 
global change and space weather. 
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